“An approach for constraining mantle viscosities through assimilation of paleo sea level data into a glacial isostatic adjustment model” by Schachtschneider et al., present a data assimilation approach to recover mantle viscosity parameters for a 3-layer glacial isostatic adjustment (GIA) model using synthetic relative sea-level change data. Based on a range of experiments, they conclude their particle filter based method was able to find the ground truth mantle viscosity parameters given different initial ensemble probabilities, observational uncertainties and spatial-temporal data coverage. This is a novel application of the data assimilation approach for GIA modelling studies, which has a good potential to help the community to better constrain the mantle viscosity parameters of a commonly-used 3-layer GIA model. This is a well-written manuscript which explains a complex subject in a very clear way. I have two major points and other minor comments. Provided the authors address these points in a subsequent revision I would certainly support the publication of this manuscript.

Major points

Firstly, given this is mainly a methodological paper for finding the mantle viscosity parameters, it is unclear what the major advantages of this data assimilation approach are compared to other commonly-used statistical approaches for finding the best-fit GIA model parameters. For example, a common approach is to use chi-square misfits to find the confidence interval of mantle viscosity (see Lambeck et al., 2014), which allows the best-fit GIA parameters to locate anywhere within the initial distribution. The authors have mentioned that their approach can converge to viscosity values that are not covered by the initial ensemble, which can also be approximately achieved given a denser sampling of possible viscosity parameter values within a forward modelling scenario. In this case, the computation time for this new approach is a very important information to report. The authors should at least demonstrate one advantage of this new approach either regarding the inversion power or the computational efficiency.

Secondly, the experiment set-up is over-idealized. In order to apply this method for real GIA problems in the future, I suggest the authors consider relaxing several assumptions, or at least discuss how relaxing these assumptions will impact the final results:

1. The authors assume a temporally-uniform relative sea level (RSL) data uncertainty of 0.5/0.25/0.1 m throughout the Last Deglaciation, which is however very difficult to achieve in reality. The RSL data uncertainty tends to be larger for earlier time intervals (e.g., a large proportion of pre-Holocene RSL data are coral-based records, which usually have >2 m uncertainty range) and gradually decreases for more recent time intervals as more stratigraphy-based data became available. Therefore, the assumption of 0.5 m upper limit for RSL data is not solid, and it is important to check how a temporally variable (sea level index point) SLIP uncertainty distribution and a larger SLIP uncertainty would impact the assimilation results.

2. The authors assume a σinit of 2x1020/2x1019 Pa s for lower/upper mantle viscosity, which is a very small range given that the authors suggest the initial ensemble should cover the ground truth parameter value (c.f. the Lambeck et al., 2014 search space is 1019 -1021 and 5x1020-1024 Pa s for upper and lower mantle respectively). My query is what if the authors use a larger σinit value, would it affect the final assimilation results?

3. The authors assume no ice loading history uncertainty in the synthetic experiments, but the ice loading history is the biggest uncertainty for GIA modelling problems. It would be useful to discuss the potential problems caused by uncertain ice loading history, in other words, are there some possible solutions for the situation if we are uncertain about where to propagate the particles?

Minor points:

-l25: The authors should pay attention to the differing role of mantle viscosity in controlling RSL change in the near-field and far-field regions. For example, in line 25 here, on-going GIA-governed far-field sea level change is primarily due to the ‘fingerprint’ effect (i.e., the spatially variable RSL change associated with changes to the shape of the gravitational field), which is largely not related to the mantle viscosity parameters. The mantle viscosity parameters are more important for near-field regions where viscoelastic land deformation dominate the local RSL change. Similarly, for line 186, the Earth deformation only dominates the near-field region, not RSL change for all regions. I would suggest the authors including some statements about the RSL change differences between the near- and far-field regions and highlight that the mantle viscosity parameter are more important for the near-field RSL change.

-l65: What is the point of this paragraph? do you use this approach or not?

-l69: Including a reference for the particle filter here would be useful.

-l93: What does the subscript i mean for equation 2?

-l130: How does the ensemble propagate through time? Do the authors run 50 forward GIA models with different mantle viscosity parameters or just one reference model to guide the propagation length and direction? I think you will update the mantle viscosity parameters during the deglaciation experiment. In order to calculate the RMS for that ensemble member, do you go back and calculate RSL value for this updated model from the start of the experiment (25.5/9.5 ka BP)?

-l146: It would be useful to mention that RSL rate is a non-standard type of observation; it is more common to compare model output with absolute RSL observations. When dealing with real-world data, the uncertainty on RSL rate will be greater than the uncertainty on absolute RSL due to the cumulative impacts of spatial-temporal uncertainty on the individual observations.  

-l167: You quote the uncertainty on RSL observations, but you compare model output with RSL rates. It would be useful to know the RSL rate uncertainty (m/ka) produced by 0.1 and 0.5 m RSL uncertainty.

-l169: Please define what you mean by ‘initial offset’ here.

-l198: Since the synthetic observations are RSL rate, why are the RMS values expressed in m instead of m/ka?

-Figure5: Can the authors slightly expand on why there are RMS spikes after meltwater pulses? Specifically, meltwater pulses are produced by large ice mass loss from North American and Fennoscandia Ice Sheets which will produce a large signal of land deformation in those regions (which could be larger than the observational uncertainty), this will help the assimilation process to find the optimal solution. Therefore, given a better converged mantle viscosity ensemble after meltwater pulses, why are there such large spikes immediately after the improvement in the mantle viscosity solution?

-l225: Both the selected regions for Setup 2 use near-field region SLIPs where RSL change are sensitive to the mantle viscosity value. If the authors just used the far-field region SLIPs, does the assimilation approach work as well? If not, the authors should note their method should only be applied to near-field GIA problems.

-l259: “As a consequence less models…”, should use “fewer models” here.

-l296: “There are no large ice mass changes after 10 ka BP.” However, based on Figure 7 of this paper, the continental ice volume was still changing after between 10 and 5 ka BP.

-l318: It would be useful if more details about computation time can be given here.

-l324: How would the temporal variability of the data availability affect the results?

-l325: “… ~8 to 12 ka BP and is considerable smaller”, should be considerably smaller here.

-l337: I suggest the authors cite Khan et al., (2019) as the reference for the HOLSEA community here.

References

Lambeck, K., Rouby, H., Purcell, A., Sun, Y. and Sambridge, M., 2014. Sea level and global ice volumes from the Last Glacial Maximum to the Holocene. Proceedings of the National Academy of Sciences, 111(43), pp.15296-15303.

Khan, N.S., Horton, B.P., Engelhart, S., Rovere, A., Vacchi, M., Ashe, E.L., Törnqvist, T.E., Dutton, A., Hijma, M.P. and Shennan, I., 2019. Inception of a global atlas of sea levels since the Last Glacial Maximum. Quaternary Science Reviews, 220, pp.359-371.

The manuscript discusses the estimation of mantel viscocities using sea-level observation in a particle filter. The manuscript is well written, and the results are interesting even though the set-up is highly simplified. I suggest publication after the following minor comments are taken into account.

Line 78: Strange sentence, not sure what the authors want to convey.

Line 79: Note that the output of the filter is the weighted posterior ensemble, from which a weighted mean can be calculated. This mean can be a poor estimate of the posterior pdf if that pdf is strongly non-Gaussian. Please add a small discussion of these facts.

Line 102: Resampling hardly changes the ensemble variance of the weighted ensemble, which is the relevant ensemble in this case. The weighting itself reduces the variance in the ensemble.

Line 105: It would be good to mention that of the three methods the second is the correct methods, and the 1st and 3rd are approximations.

Line 111: I assume the authors mean N(0,a \sigma^U,L). Note that it is common to have the (co)variance of the distribution as the second argument of N(..,..) and not the standard deviation.

Line 114: Setting a=0.5 \sigma is a large value to perturb each viscocity with. However, it is perhaps good to mention that the standard deviation of the ensemble as a whole only increases by a factor 1.12 (= sqrt(1 + 0.25)) through this procedure.

Line 115: Probabilistic resampling means that one has to draw N random numbers, where N is the ensemble size. A more accurate resampling method is Stochastic Universal Resampling, in which only one random number is drawn, and it is also faster! This could be mentioned. Since after resampling relatively large random perturbations are added to the particles the difference will be minor in this case.

Line 165: mu

Line 174-175: Please remove this sentence, it is just a repetition of what was said before.

Since the variance in the initial ensemble is chosen as large as the mean many initial particles will have negative viscocities. What is done when a negative viscocity is drawn? A similar question for later in the run, what is done if the 'jittering' after resampling produces negative viscocities? Related to this, the figures show different mean viscocities than the table 2. Please correct the one that is incorrect.

It would be good to show the effective ensemble size, defined as N_eff = 1/sum_i (w_i^2) in which the w_i are the normalized weights, such that sum_i w_i = 1. This allows the readers to judge the quality of the ensemble. My suspicion is that N_eff is rather low, as low as 2-5 members at time, which is close to degeneracy.

Fig 10: caption, change left and right to top and bottom.

Line 309: I don't understand this sentence, please clarify.

Section 2: I'm not an expert in this field and would suggest providing the set of equations being solved to gain an idea of the complexity of the problem at hand.

Section 7: It would be good to also include a discussion of the accuracy of the underlying ice model and its expected influence on the results.

The authors use a particle filter for parameter estimation in a visco-elastic model of the lithosphere and mantle. The goal here is to find a set of parameters in a model that leads to a solution consistent with a set of observations. In this work the unknown parameters are the viscosities of the lower and upper mantle. The authors study the effect of different parameter initializations and observation uncertainty on the performance of the filter with clearly stated results and conclusions.

The paper is well organized and well-presented. I recommend publication but suggest that the authors address the following niggles:    

Line 78 I don’t get this sentence “the ensemble members did not mix they were called particles”

Line 85 This is also a very good review paper:

PJ Van Leeuwen, HR Künsch, L Nerger, R Potthast, S Reich, Particle filters for high‐dimensional geoscience applications: A review, Quarterly Journal of the Royal Meteorological Society 145 (723), 2335-2365, 2019

To find more generic references on particle filters, I suggest that the authors look at: 

Doucet, A.; Johansen, A. M. A tutorial on particle filtering and smoothing: fifteen years later. The Oxford handbook of nonlinear filtering, 656–704, Oxford Univ. Press, Oxford, 2011.

Line 110 An alternative method to update the parameter set is to use jittering. For details:

Crisan, Dan; Míguez, Joaquín Nested particle filters for online parameter estimation in discrete-time state-space Markov models. Bernoulli 24 (2018), no. 4A, 3039–3086.

Line 138 In section 4.2 It would be good to explain what the observations are, the choice of the distribution for observation noise (the observation uncertainty parameter is stated at line 165 but not the distribution) ,  and give the  expression for the parameter likelihood function.

Finally, perhaps you can put the VILMA equations in an appendix.

This manuscript introduces a combination of a data assimilation method with a GIA model. The new combination approach is supposed to, ideally, better determine mantle viscosities or, in future, 3D earth parameters. To achieve this, a set of global synthetical relative sea level rates is generated. The authors perform a number of tests to convince the reader about the favourable outcome of their new approach. They seem to be new in the field of GIA modelling. I am not aware of any previous GIA work with the exception of Volker Klemann, who has a strong background in this field.

The study is very interesting and the approach may receive much interest in the GIA and earth rheology communities. The text is well written and a smooth read. Figures and tables are clear and support the text.

Nonetheless, I am somewhat disappointed about the whole manuscript. Main reasons are that, despite the nice presentation and different tests made by the authors, (1) nothing is presented about the performance of the new approach compared to 'GIA standard' methods, (2) no interesting conclusions are drawn that present a step forward in GIA/earth rheology research, and (3) it is a very idealized experiment because the synthetic data does not represent typical data used in GIA modelling. RSL rates is not used, moreover, rates are hard to find in real data as there is rarely a large number of samples at a certain location available. Also, uncertainties of rates are much larger than used in these setups. Overall, the main message is that the technical combination of two codes is working and gives results that are expected. The results thus, at this stage, cannot help to further advance our understanding of GIA or earth parameters.

However, I think the manuscript can be elevated if my main concerns can be addressed.

(1) Especially, I would like to see a comparison with a 'standard' GIA investigation, where modelled RSL rates from a set of pre-defined models (e.g., 50 models covering the viscosity ranges in your experiment) is compared to the synthetic set and the misfit is determined so that a best model of such set is identified. Is the best-fitting model comparable to the final assimilation model? Which misfit is better? What is the computation time for both approaches? At which point is it better to use the assimilation approach? This would help the reader to get more perspective if this approach can help advance our understanding of GIA and the determination of earth parameters.

(2) I miss new findings or hints that can help the community. The manuscript presents the approach with some tests, which gives it the style of a technical note rather than a scientific study. The authors should at least present 1 or 2 major conclusions that can be drawn from the tests.

(3) The reliability of the rates set should be further discussed in comparison to real world data. You mention some shortcomings but they are not put into perspective with real data availability. How many locations can actually give you solid RSL rates? What is a realistic error of such RSL rates? This should definitely be addressed as RSL data are concerned with time and height errors. You did not include time errors which are actually much larger! Are there enough locations with rates at times where there is a strong RSL fall? Such discussion would help the reader to get more insight on the reliability and evaluate the success of your approach.

(4) A discussion is needed on the tested parameters. Just analyzing two mantle viscosities is very idealized. There is a trade-off between the thickness of the lithosphere and mantle viscosity. The reader should be informed. Similarly, a note on ice model uncertainty and its potential impact on the results should be added.

Minor remarks

The paper is written from a quite technical perspective. In the introduction, focus is a lot on the assimilation approach but I would like to see a paragraph from the 'GIA site' with an overview of previous attempts to get more insights from GIA modelling with alternate approaches. Studies by Steffen & Kaufmann (2005), Al-Attar & Tromp (2013) and Caron et al. (2017) should help here. Similarly, the discussion does not contain much references to other works. Are all these findings/conclusions new?

References

Al-Attar, D., Tromp, J., 2013. Sensitivity kernels for viscoelastic loading based on adjoint methods. GJI, doi:10.1093/gji/ggt395.

Caron, L. et al., 2017. Inverting Glacial Isostatic Adjustment signal using Bayesian framework and two linearly relaxing rheologies. GJI, doi: 10.1093/gji/ggx083.

Steffen H., Kaufmann, G., 2005. Glacial isostatic adjustment of Scandinavia and northwestern Europe and the radial viscosity structure of the Earth’s mantle. GJI, doi:10.1111/j.1365-246X.2005.02740.x.