07 Jun 2021
07 Jun 2021
An approach for constraining mantle viscosities through assimilation of paleo sea level data into a glacial isostatic adjustment model
 ^{1}Helmholtz Centre Potsdam GFZ German Research Centre for Geosciences, Telegrafenberg, 14473 Potsdam, Germany
 ^{2}Freie Universität Berlin, Kaiserswerther Str. 1618, 14195 Berlin, Germany
 ^{1}Helmholtz Centre Potsdam GFZ German Research Centre for Geosciences, Telegrafenberg, 14473 Potsdam, Germany
 ^{2}Freie Universität Berlin, Kaiserswerther Str. 1618, 14195 Berlin, Germany
Abstract. Glacial isostatic adjustment is largely governed by rheological properties of the Earth's mantle. Large mass redistributions in the oceancryosphere system and the subsequent response of the viscoelastic Earth have led to dramatic sea level changes in the past. This process is ongoing and in order to understand and predict current and future sea level changes the knowledge of mantle properties such as viscosity is essential. In this study we present a method to obtain estimates of mantle viscosities by assimilation of relative sea level data into a viscoelastic model of the lithosphere and mantle. We set up a particle filter with probabilistic resampling. In an identical twin experiment we show that mantle viscosities can be recovered in a glacial isostatic adjustment model of a simple three layer earth structure consisting of an elastic lithosphere and two mantle layers of different viscosity. In two scenarios we investigate the dependence of the ensemblebehavior on the ensemble initialization and observation uncertainties and show that the recovery is successful if the target parameter values are properly sampled by the initial ensemble probability distribution. This even includes cases in which the target viscosity values are located far in the tail of the initial ensemble probability distribution. We then successfully apply the method to two special cases that are relevant for the assimilation of real observations: 1) using observations taken from a single region only, here Laurentide and Fennoscandia, respectively, and 2) using only observations from the last 10 kyrs.
Reyko Schachtschneider et al.
Status: final response (author comments only)

RC1: 'Comment on npg202122', Anonymous Referee #1, 25 Jun 2021
“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 3layer glacial isostatic adjustment (GIA) model using synthetic relative sealevel 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 spatialtemporal 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 commonlyused 3layer GIA model. This is a wellwritten 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 commonlyused statistical approaches for finding the bestfit GIA model parameters. For example, a common approach is to use chisquare misfits to find the confidence interval of mantle viscosity (see Lambeck et al., 2014), which allows the bestfit 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 setup is overidealized. 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 temporallyuniform 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 preHolocene RSL data are coralbased records, which usually have >2 m uncertainty range) and gradually decreases for more recent time intervals as more stratigraphybased 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 5x10201024 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 nearfield and farfield regions. For example, in line 25 here, ongoing GIAgoverned farfield 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 nearfield regions where viscoelastic land deformation dominate the local RSL change. Similarly, for line 186, the Earth deformation only dominates the nearfield region, not RSL change for all regions. I would suggest the authors including some statements about the RSL change differences between the near and farfield regions and highlight that the mantle viscosity parameter are more important for the nearfield 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 nonstandard type of observation; it is more common to compare model output with absolute RSL observations. When dealing with realworld data, the uncertainty on RSL rate will be greater than the uncertainty on absolute RSL due to the cumulative impacts of spatialtemporal 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 nearfield region SLIPs where RSL change are sensitive to the mantle viscosity value. If the authors just used the farfield region SLIPs, does the assimilation approach work as well? If not, the authors should note their method should only be applied to nearfield 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.1529615303.
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.359371.
 AC1: 'Reply on RC1', Reyko Schachtschneider, 04 Oct 2021

RC2: 'Comment on npg202122', Peter Jan van Leeuwen, 27 Jun 2021
The manuscript discusses the estimation of mantel viscocities using sealevel observation in a particle filter. The manuscript is well written, and the results are interesting even though the setup 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 nonGaussian. 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 174175: 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 25 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.
 AC2: 'Reply on RC2', Reyko Schachtschneider, 04 Oct 2021

RC3: 'Comment on npg202122', Dan Crisan, 15 Jul 2021
The authors use a particle filter for parameter estimation in a viscoelastic 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 wellpresented. 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), 23352365, 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 discretetime statespace 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.
 AC3: 'Reply on RC3', Reyko Schachtschneider, 04 Oct 2021

RC4: 'Comment on npg202122', Anonymous Referee #4, 11 Aug 2021
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 predefined 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 bestfitting 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 tradeoff 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), AlAttar & 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
AlAttar, 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.1365246X.2005.02740.x.
AC4: 'Reply on RC4', Reyko Schachtschneider, 04 Oct 2021
Dear reviewer,
thank you for your constructive comments and questions about our manuscript. In the attached pdf file you can find our detailed answers to your comments.
Kind regards
On behalf of the authors, Reyko Schachtschneider
 AC5: 'Reply on AC4', Reyko Schachtschneider, 04 Oct 2021

AC4: 'Reply on RC4', Reyko Schachtschneider, 04 Oct 2021
Reyko Schachtschneider et al.
Reyko Schachtschneider et al.
Viewed
HTML  XML  Total  BibTeX  EndNote  

663  103  25  791  5  6 
 HTML: 663
 PDF: 103
 XML: 25
 Total: 791
 BibTeX: 5
 EndNote: 6
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1