see attached pdf

general comments

The manuscript “Simulation of a fully coupled 3D GIA - ice-sheet model for the Antarctic Ice Sheet over a glacial cycle” by Caroline J. Van Calcar et al. model the influence of GIA feedback on evolution of the Antarctic Ice Sheet over the full glacial cycle from last interglacial to present. This study uses state-of-the-art approaches and models such as incorporating the lateral heterogeneity in the mantle structure using a recent seismic model of Lloyd et al. 2020 to develop two different 3D Earth structures assuming a composite rheology and applying a dynamic coupled GIA-ice sheet model with variable coupling time steps. Therefore, the results and methodology will have a significant and excellent contribution to the community. However, it needs more clarification and major revisions to be considered for publication. 

- Some aspects of the applied methodology are not clear and need further explanation. One important aspect is the GIA feedback and the results of the GIA model that is used for the coupling approach. Throughout the manuscript the authors use the term “deformation” for such model output, while the GIA signal is a combination of gravitational and rotational effects as well. Therefore, the surface deformation is not the only factor that affects the ice evolution, but also the changes in gravity field (geoid) and consequently the sea level which essentially defines the topography. This needs to be clarified if (and why) the applied method has ignored such effects.

- The approach to infer the viscosity field is still not clear (Line 230-290). A detailed description is required on how the absolute viscosity values (the term stated in Line 281) are inferred from seismic anomalies.   

- The elastic thickness of the lithosphere considered for the 1D models differ from that of the 3D models, the former is set to 100 km and the later is set to 35 km (Line 285). Please explain the reason for choosing these two different values. The lithosphere thickness is an important parameter particularly where there is localized loading and deformation, as a thicker lithosphere can have a damping effect on the GIA signal. So, to have a fair comparison, this parameter should have an equal value (or very close values) in both 1D and 3D settings. Therefore, at least one additional set of simulations is required where the elastic thickness of the lithosphere is set to 35 km in the 1D models. The associated model results should be revised based on this simulation. Furthermore, as mentioned in the previous comments, 35 km is not a good representative of the lithosphere thickness across East Antarctica. So, a sensitivity analysis with a thicker lithosphere (e.g., 100 km as set in the 1D models) in the 3-D simulations would be beneficial.

- The model outputs that are presented in the result section are based on only one iteration over the coupling time step, as the maximum differences are flagged as the outliers, and it is declared that the absolute mean of the maximum differences is 2.4 m. However, based on Fig. 5, there seems to be a significant difference between changes in the ice thickness and bed elevation of the second iteration and those of the first iteration. Also, the comparisons that are provided in section 2.4.3 are not clear that belong to which simulations. Therefore, this section requires more explanation and clarification with supporting figures in the supplementary section. Furthermore, based on the provided reasoning and Fig. 5, it is not convincing that only 1 iteration is sufficient to achieve the convergence criteria. In addition, if the results are based on only one iteration, section 2.4 (including Fig. 4) and elsewhere in the manuscript needs to be revised to declare that only one iteration per coupling time step is performed. Otherwise, it would be misleading for the reader.  

- The result section can be developed further to include detailed discussion about different aspects of the model results, particularly when compared with other studies that performed GIA-ice sheet coupling approach. There is lack of model comparison with observational constraints, some constraints that can be considered here include: present-day (grounded) ice elevation and grounding line position, surface exposure age data, and relative sea-level records. As mentioned, while the similarities between the model results of this study and other coupled GIA-ice sheet studies are considered, there is lack of investigation of the mismatches, data-model comparison enables the authors for performing a better evaluation in this regard. Furthermore, the current plots do not adequately support some of the provided information/numbers within the text.

- As the final remark, the dynamic coupling is introduced as a feature element of the modelling in this study in comparison to other studies with coupled simulations (e.g., Line 337-340). However, its implications are not explored thoroughly. It would be interesting to see how the results differ from the case when coupled GIA-ice sheet model is done alternately at once over the entire simulation time in order to highlight its impact.  

specific comments

Line 17 in abstract. “ … on a high temporal resolution …”, since a combination of temporal resolution is used in the manuscript from 5000-500 years I would suggest not using the term “high temporal resolution” and please provide an accurate statement accordingly.      

Line 19. Please be specific about the applied seismic model.

Line 21. Please be specific about the region(s) where these maximum differences occur.  

Line 26 in introduction. “ …. is the response of the solid Earth  …”, one important aspect of GIA is its effect on the gravity field of the Earth which is essentially tied to the definition of sea level (and so the bedrock elevation). Therefore, solid Earth does not reflect this point adequately. Please revise this sentence accordingly.   

Line 34-36. “ … causing a local shoaling of water and an outward movement of the grounding line to position p3 (Fig. 1). As a consequence, the GIA feedback slows down migration of the grounding line (…) and acts as a negative feedback (e.g. Konrad et al., 2015)”. The relation between ice flux, water depth and the ice thickness need to be clearly stated here as the reason behind the stabilizing effect of the GIA and the readvance of the grounding line. Also, is there a reason that the final sentence and its associated reference is not followed immediately after the previous sentence with a combined list of references?  

Figure 1. Change Elastic “crust” to Elastic “lithosphere” as the latter also includes the part of the upper mantle that is brittle in addition to the crust.

Line 45-48. The gravitational effect of the GIA impact the grounding line migration and the ice retreat, the phrase mentioned “apart from the effect on the grounding line” is confusing and does not apply in this case. Please revise accordingly.

Line 50. Please clarify how the final feedback “Finally, ..” differ from the process illustrated in Figure 1 and the beginning of the introduction section?

Line 65. “Some 1D GIA-sea level models also account for relative sea level change”, please clarify how these differ from the 1D models that account for gravity field perturbations and displacements as mentioned in previous sentence.

Line 65-66. The models and references cited here are referring to the studies rather than different models. For example Nield et al., essentially uses Whitehouse et al., 2012 GIA model and similarly DeConto et al. and Gomez et al. uses a similar GIA model. So, this is confusing whether the authors are using examples for different studies or different computational methods. If latter, there are some references that can be referred to as the first generation of the 1D GIA models that should be used. Please clarify.

Line 68. The relevance of the reference “Geruo et al., 2013” is not clear here, please clarify the reason for choosing this reference.

Line 70. “although for the Eurasian ice sheet, … (van den Berg et al. 2008)”, the relevance of this sentence to previous sentence is not clear and also the cited reference does not include such conclusion. Please clarify and also mention which page and paragraph this statement is mentioned in van den Berg et al. 2008.

Line 103-104. “from the previous interglacial”, please indicate clearly the temporal coverage with respect to the present time.

Line 109. “In this study we neglect the spatial variations in sea level.”, how is sea level defined here? This makes a confusion as the GIA model is essentially referred to as a “sea-level solver”. Please clarify this and add more description.

Line 122 in method. “individually or simultaneously on different equidistant grids for each ice sheet”. How does simulation on different grids differ from the term individually mentioned here? Please rephrase this sentence for more clarity.

Line 123. Does the usage of the term “typically” mean that the horizontal resolution can be adjusted, e.g., to higher resolution over a given region? If there is an option for adjusting the grid resolution, what is the justification for choosing a 40-km grid resolution over AIS and not 20 km? Are there any technical difficulties involved? Please clarify.

Line 136. Indicate here that the model does not take into account this regional sea-level variations.  

Line 136-137. By stating “similar throughout Antarctica”, does this mean the effect of the NH ice sheet on the evolution of Antarctica? If so, this effect has been recognized to be noticeable by the study “Antarctic ice dynamics amplified by Northern Hemisphere sea-level forcing”, Gomez et al., Nature, 2020. Please clarify.

Line 137. “The effect of the AIS itself on regional sea level is more important”, what aspect is considered here when stating “more important”, as the change in regional sea level provides feedback to AIS and this is the foundation of the coupling approach which is considered to be quite important.   

Line 140. “the effect of regional sea level variations is a second order effect compared to the GMSL variations …..”. What is the justification for this statement?

Line 136-141. “The effect of the northern … yet included in this model”. This part seems to include information that are not necessary and can be removed.

Line 147. “at time steps of 1 year”, is 1 year the computational time step of ANICE? Please mention this temporal resolution of the ice sheet model at the beginning of this section where the spatial resolution of 20 km and 40 km are noted.

Line 159. “it computes bedrock …. on a spherical Earth”,  what are the assumptions of the Earth model here? Does the model assume it as an incompressible material? If so, please add this information. Also, what is the assumption for the rheology of the Earth?  

Line 11. “over the full glacial cycle”, indicate the time window of this study. 

Line 174. Wu et al. (2004) demonstrate that achieving the solution to include self gravity below 1 per cent error requires 4-5 iterations. Please quantify how much the lower number of iterations affect the accuracy of the results.   

Line 175. “The same iteration within each time step …”, what is the purpose of this, please clarify and also state that whether this has been performed in this study or not.

Line 190-194. Add relevant references and relative plots/results for justification of the statements.  

Line 196-197. Is the lithosphere thickness considered as a constant value? Please provide information about the applied lithosphere thickness value/model.

Line 201. “Since the difference in deformation is insignificant, …”, the maximum difference of the performed test is associated with a viscosity of 10²⁰ Pas. This is while there are regions of very low viscosity beneath western AIS. It would be useful to see the similar graph for a lower viscosity of e.g., 10¹⁸ Pas. Also, the information about the radius of the applied load in the test is not provided in the description of Fig. S2, please add this information.

Line 218. Cite table 1.

Line 218. The order of the details is making confusion. I would suggest to bring the details according to the depth, as in table 1, starting from the lithosphere, so one would not need to wonder about the viscosity values considered for the upper parts of the mantle above 210 km of depth.

Line 220. Whie considering a thin lithosphere is appropriate for the western section of AIS, the eastern region is characterised by thicker lithosphere. A sensitivity test is required to explore the effect of considering a thick lithosphere on the results

Line 221. “..with specific dislocation and diffusion creep parameters..”, the assumption of a composite rheology should be noted in section 2.2 (see comment of Line 159).

Line 226. Table 1, change crust to lithosphere (as crust is a compositional layer of the Earth).

Line 251. If A is a constant, it should be mentioned here.

Line 254-255. Which seismic models? move this information that is mentioned later to this sentence. Provide details about approach 3 in Ivins et al., 2021, if provided it is not clearly linked in the following sentences.

Line 272-274. It would be useful to have a figure in the supplementary material showing the viscosity profile for the two models, 3Ddry and 3Dwet, at some locations in e.g., Amundsen Sea sector, Antarctic Peninsula, etc.

Line 282. “…assuming only a temperature profile and not a viscosity profile.”, it is not clear what assumption is made here.

Line 230-290. Please cite Table 1 where there are numbers pointing to the values in this table.  

Line 289-290. “Therefore, the rheology …”, this sentence is implicit in the previous sentence so it can be removed.

Line 301-302. What are considered as the ANICE output here? The topography and the ice loading? Please specify.  

Line 310. “Furthermore, the ice thickness is linearly interpolated …”, this is confusing which method is used for interpolation of the ANICE output on the GIA grid, as in Line 306 it is stated that the quadrant method is used for gridding from a coarser ANICE grid to a somewhat finer grid of the GIA FE model. Please clarify.

Line 298-311. This section (section 2.3) can be placed in the supplementary material and a sentence can be added in the coupling section to mention that the GIA and the ice model outputs are generated on different grids and the corresponding interpolation method is described in the supplementary section.  

Line 312-345. In my opinion, instead of explaining the coupling approach for each coupling time step, first the big picture of the modelling can be described, including the coupling scheme. This is because the current explanations of the model make some confusion that will be resolved later but would be better if the steps are clear as the reader goes through the manuscript. For example, in GIA modelling, the difference between the observed modern topography and the predicted one is used to update the initial topography, and this is done in an iterative process. This is mentioned later (lines 342-343), but it would be much clearer if the entire modelling scheme is shown in the figure and included in the bullet points.

Line 335. Figure 4. Following my previous comment, as there are n steps involved, the figure cuts just after coupling time step 2. I would recommend revising the figure so that it shows the entire performed approach including the following time steps (or the final, n^th time step) after time step 2 and the iteration for the convergence of the present-day topography.  

Line 315-316. And also in Lines 145 and 187. “ …. developed that alternates between the models per time step of 500 to 5000 years.”, from the statements it is not clear that the model alternates between these coupling time steps, it may give an impression that the model is performed for each of these coupling time step. So, please revise the associate sentences accordingly. In addition, please provide a short sentence of the criteria used for choosing the time step.

Line 349. Fig 5 can be moved to the supplementary section.

Line 348-349. For that given rheology? Is the number of iterations also dependent on the assumed rheology (e.g., regions with lower viscosity experience larger GIA signal upon ice load change)?  

Line 350. Mention the associated Earth model (e.g., 1D21 rheology).  

Line 347-350. I would recommend starting the paragraph with a general statement, e.g., the number of iterations per coupling time step to converge is dependant on the rheology and the size of the ice load. Our simulation shows that for a given rheology, the 1D21 model, The coupled model requires three iterations per coupling …. or something like this.

Line 365-366. “the uncertainty range of the GIA FE model based on uncertainties from the rheological model such as background temperature and seismic velocity …”, using “such as” here do not seem a proper link here. I would suggest the usage of other terms such as associated with, etc.   

Line 379-385. These sentences can be reordered for the purpose of clarity. Here is my recommendation: First bring up the point that “The convergence of the coupled model depends on the length of the coupling time step, since smaller time steps increase the number of grid cells converging to zero.” Then mention that there are other factors than the length of the time step that affect the convergence: “The convergence is also highly dependent on the change in deformation and ice thickness such that the time steps need to be chosen sufficiently small to have nearly linear changes in ice thickness and bedrock elevation.” Then declare the disadvantage of smaller time steps: “On the other hand, a large time step is desirable to limit …”

Line 389-390. “but their method assumes a constant topography during one coupling timestep which requires smaller timesteps than the coupling method presented in this study.”, what does it mean by “a constant topography” in the study of Han et al., 2022? As stated in section 2.4, a total deformation from GIA model is passed to the ice sheet and the next coupling time step, which seems a similar concept as the coupled simulation applied by Han et al. Please clarify. Also, Han et al. (2022) suggested 200 years as the shortest and preferred coupling time interval for glacial-cycle simulations. Please revise the sentence. Considering the 0.2 kyr as their preferred coupling time interval, how do you justify the selection of the 500 year coupling time interval?  

Line 398-401. “The absolute maximum …” it is difficult to follow these sentences, Are the two numbers associated with the absolute maximum difference in ice thickness between the 1-iteration and multi-iteration simulations? Please clarify. “the absolute mean of the maximum differences …. Is 2.4 m”, To which parameter this number belong?, “The maximum difference in ice thickness at present-day …”, how much is the difference and where it occurs?  

Line 414. The term “subsidence” would be a more appropriate choice here.

Line 414. “Ice shelfs in West Antarctica will melt less … due to the higher bedrock elevation”, any reference?

Line 416-418. “Differences in ice sheet evolution during the deglaciation phase are then mainly caused by a different topography at last glacial maximum rather than differences in rheology.”, It is difficult to follow why the previous sentence leads to this conclusion. Also, the relevance of this statement to the following sentences are not clear. Please clarify.

Line 432-500. The results section can be presented in another way, by categorizing the subsections based on different time intervals. For example, one subsection is assigned to LIG to LGM, another subsection covers the last deglacial retreat, either from LGM to present day or this can be divided into LGM to late Holocene, and late Holocene to present day. This approach facilitates investigating the model outputs.

Line 444. In “At present day, the ice is up to 1 km thinner around the grounding ….”, this information cannot be inferred from the figures provided, so it is suggested to provide a zoomed map of the places where such information are provided. Also, a figure can be added as “data-model” comparison where the present-day observed ice thickness and grounding line position is compared to the model results.

Line 438. “These simulations also allow to study the differences between ..”, this sentence seems that do not follow the flow of the text and can be removed.

Line 451-501. As mentioned in earlier comments of the result section, the model outputs provided in this subsection can be compared to observational constraints, such as the modern ice thickness, to better investigate the efficiency of the developed 3D models.

Line 454-454. Any statement should be cited to the associated figure, if not provided, one should be generated to support the statement within the text.

Line 454. “In the 1D simulations, the bedrock subsides approximately 500 meter less”, the same comment as the previous one, the statements should be supported by appropriate plots.

Line 458-459. “During the deglaciation phase, the Ross and Filchner-Ronne Ice Shelfs retreat fast due climate warming, similar to other studies of the AIS evolution suggest (e.g. Albrecht et al., 2020).”, As mentioned earlier in the general comments of the result section, beside the similarities, the differences in the model results should also be investigated. As an example, there are some discrepancies on the timing of the thinning around ASE and Ross embayment, some studies suggest earlier deglaciation while in the study referred here the major deglaciation occur after 10 ka (though the more appropriate referencing is the part 2 of Albrecht et al., 2020, https://tc.copernicus.org/articles/14/633/2020/). So a more detailed comparison and discussion is required.

Line 464. “Using a 3D viscosity leads to a difference in grounding line position of up to 500 km and a difference in ice thickness of up to 1.5 km at present-day”, this statement cannot be inferred from Fig. 7. Revise the plot accordingly. Also, “up to 500 km”, in which region?

Line 470. Figure 7. the comparison should include all the models involved, particularly the two different estimated 3D rheologies, i.e., also add the 3Dwet results: 1D20-3Dwet and 1D21-3Dwet.

Line 473-474. “the maximum difference in grounding line position is approximately 40 km (Fig. 8.g)”, in which region? This statement cannot be inferred from the plot.

Line 490. Figure 9. change the plot based on the unit km³.

Technical corrections

Line 14 in abstract. “Most studies assume a relatively high laterally homogenous response time of the bedrock”. I would recommend using other adjectives instead of “high” that are used in reference to the “response time”, such as slow or long.

Line 17. Please revise the sentence “The feedback effect into account …” as it seems that the word “takes” is missed.  

Line 18. Please change FE to finite element  

Line 21. Change “…, to differences in ice thickness …” to ““…, and differences in ice thickness …”

Line 31 in introduction. Put “,” before respectively.  

Line 42. “the dashed brown the new bedrock surface”, the word “is” is missing.

Line 59. Change Earths to Earth or Earth’s

Line 62. “Another approach …”, the verb should be singular. Revise accordingly, e.g., Another approach to compute GIA is using …

Line 82. “The only model that coupled 3D”, a past tense is used here while in some other places the present tense has been used. Please be consistent time wise when referring to a study.

Line 106. If FE stands for finite element, it has not been introduced before. Please indicate the full name of the method where first appeared in the text.   

 Line 112. I would recommend rephrasing of “as 3D Earth structures” as: “in comparison to 3D Earth structures“.

Line 114. Chane to “This method … ”

Line 140 in method. “sea-level” variations

Line 161. Can be written as: “…. based on ABAQUS is its flexibility as its grid size…”

Line 162. “and FE models operate in the time domain …” I would suggest to make this part as a separate sentence, otherwise the whole sentence is long and more difficult to follow.  

Line 219. Change to consistent

Line 252. In Tx,y, the x,y indices should appear the same as the equation (as subscripts)

Line 279-280. “who obtained viscosity by scaling seismic …” and “background viscosity can be obtained from …”, using the term “obtain” in these sentences give the impression that the viscosity values are determined or tuned rather than being modelled. I would recommend using alternative terms such as model, estimate, infer, etc.

Line 284. “two experiments …”, I would recommend to rephrase this sentence for clarification, something like: “Two experiments are performed using a 1D rheology with two different upper mantle viscosity profiles: 10²⁰ Pa.s (hereafter referred to as 1D20) and 10²¹ Pa.s (hereafter referred to as 1D21). These values are consistent with the lower and upper boundaries… . The elastic lithospheric thickness is the same for both experiments and is set to 100 km. “. A general comment: please keep the sentences short so that it will be easier for the readers to follow.

Line 359. It may read better:” In this case, both ice thickness and deformation at these … “

Line 364. The coupled model “converges” within an acceptable ….

Line 365. No need for comma after the GIA FE model.

Line 365-367. I recommend to move “accuracy of paleo sea level records” before “the uncertainty range of the GIA FE model” for clarity. So, the sentence reads: “This is within the accuracy of paleo sea level records and the uncertainty range of the GIA FE model …”.

Line 369-371. These two sentences can be merged as one single sentence, something like: To decrease this uncertainty, the average deformation of the last two iterations is used as the final deformation to simulate ANICE for the final iteration of a given time step.

Line 386. “ice sheet is slowly increasing till LGM,…”, increasing in what aspect? Extent? Volume? Please specify.

Line 398. Also everywhere else in the manuscript: change meter to “m”

Line 399. And elsewhere in the manuscript. Change kyears to kyr, kyr before present can be written as ka for simplicity   

Line 420. “If they are not, it is assumed that initial topography is in error.”, this is implicit in the previous sentence and can be removed.

Line 426. And also Line 429. It seems that there is an extra parentheses in the subscript of H_b,ALBMAP)

Line 434. Change to “extent”

Lin2 449. e-“h”?

Line 460. Grounding line

Line 470. Figure 7 (and also Fig. 8). Indicate the depth associated with the 3-D model viscosity map.