Review of « A microstructure-based parameterization of the effective, anisotropic elasticity tensor of snow, firn, and bubbly ice » by Sundu et al. in The Cryosphere

Summary : 

The effective stiffness tensor of snow, firn, and bubbly ice is controlled by the density, morphology, and elastic properties of the ice matrix. This control was previously studied and parameterized independently for different ranges of density: for snow (rho in [30, 500] kg/m3), firn (rho in [500, ~800] kg/m3), and bubbly ice (rho in [~800, 915] kg/m3). Here, the authors developed a new parameterization of this control that is valid on the full density range. They use the formal anisotropic Hashin-Shtrikman upper bound as a predictor of the stiffness tensor in empirical fit based on 395 finite-element simulations on tomographic images.

Main comments :

This article constitutes a valuable contribution to The Cryosphere with a sound methodology and interesting results. The paper is fairly well written to follow the work (typos to be corrected through proof review). Even if the new parameterization does not substantially outperform existing parameterizations in their porosity range of validity, it is valid on the full density range from snow to bubbly ice and does not exhibit artificial and arbitrary transition zones with density. In particular, the existing parameterizations for snow predict an effective stiffness larger than the one of ice if applied on high-density samples. Besides, the anisotropy is directly captured by the Eshelby tensor, which does not require additional fitting when using the 2-parameters fit (same beta and eta for all components) for the whole tensor. This new parameterization comes at the cost of a more complex implementation, particularly of the upper bound C_U using the Eshelby tensor derived from the correlation length of the structure. However, its computation expense remains far lower than a full finite-element simulation on the snow microstructure. 

However, I have certain comments that would need to be addressed before publication :

1. One goal of the presented work is to provide a new parameterization of the elastic tensor valid from snow to ice. To use it, one must compute the density, the correlation lengths of the given sample microstructure, the associated 4th-order Eshelby tensor, the corresponding Hashin-Shtrikman upper bound, and eventually, the empirical fit, and to juggle between Voigt and tensorial notations. The authors should provide the functions (e.g., Python or Matlab style) so the community can easily re-use this fit. Otherwise, I fear that simple density parameterization will remain the norm. The shown material must be enough to re-implement the fit, but it is prone to errors and headaches.
2. The computation of the effective isotropic transverse elastic tensor from finite-element simulations is not described in enough detail. 
1. First, what sample size (mm) and boundary conditions are used? Indeed, the convergence of the apparent sample properties into effective material properties with the simulation volume depends on the applied boundary conditions and sample density. In particular, the low-density samples of Alp-DIV likely deviate from the proposed parameterization because of the too-small sample size (Fig. 2). With this information, the robustness of the simulations can be evaluated.
2. The isotropic transverse tensor is estimated from 5 load states (Sec. 3.4) by finding the five independent components of C that minimize the L2-norm of sigma-C:epsilon. The five load states are not described. It is unclear whether a bad choice of these load states may favor better approximations of certain components when approximating the full tensor to the isotropic transverse one. What is the difference between getting the full tensor (21 components) based on 6 unit load cases and taking the theoretically non-zero components under the assumptions of transverse isotropy (e.g., Wautier et al., 2015)? In addition, the assumption of transverse isotropy makes sense for snow (deformation by gravity generally aligned with temperature gradient), but is it relevant for bubbly ice on ice sheets that may also flow in a certain horizontal direction?
3. The different models were fitted on the simulation data using a log-transformation of the elastic tensor component with a least squares regression. The density distribution of the samples is not uniform in the full density range. In particular, around 80% (?) of the samples exhibit a density between 250 and 500 kg/m3 (Tab. 1, Fig. 2). Besides, some data are highly correlated because they belong to the same time-series. The collected is already huge and the largest so far to my knowledge; however, could you discuss this point? Can we rely on this parameterization for any collected snow data, or is the fit impacted by the sampling?  Moreover, the improvements of the new parameterization do not show up in the regression coefficient (Tab. 2) or the scatter around the predictor (phi or C^U in Fig. 2).   Sampling.
4. The authors state that « the limit of φ → 1 the microstructure must tend to an isotropic state » (l.160-161). I disagree with the statement or I have not understood it. Bubbles in ice may be very flat and tend to, for instance, horizontal micro-cracks (porosity tends to zero, but anisotropy can remain constantly high). This point motivated the choice of the HS bound as a predictor but there is no prior reason for that. It appears that the collected samples (Fig. 6) of high density (phi > 0.7) are also characterized, but the sampling may be too limited to draw definite conclusions on the structure anisotropy at high density. Moreover, Fig. 6 is based on this specific feature of the HS bounds. It shows that the anisotropy of the bubbles does not affect the anisotropy of the elastic tensor. I am not convinced this is sound. Please clarify.
5. The elastic tensor depends on density as a power law with an exponent in [3, 5]. An error of 5% on density may cause an error of 15% to 25% on the elasticity components. Measuring density, even with tomography, is subjected to errors in this order of magnitude  (e.g., Proksch et al., 2015; Hagenmuller et al., 2016). The « relative » error due to anisotropy should be discussed with respect to the errors on density and not shown as the main source of uncertainty. 

Minor comments :

l11 : « the crystallographic anisotropy » ->  « to the maximal theoretical crystallographic anisotropy .» Indeed, your estimation of crystallographic anisotropy is very rough.

l22 : « the last example … » -> « Schlegel et al. have stressed ».

l24-26 : « ice matrix geometry … crystallographic orientation ». There are references for geometrical anisotropy but no for crystalline anisotropy.

l26-28: « fabric is low/high ». What does it mean? Anisotropy is high /low?

l29: « recent work wave propagations » -> ? « Hellmann et al. (2021) measured wave propagation on glacier ice and suggested … »

Figure 1: The range of density on which the existing parameterizations are supposed to work (according to their respective authors) is never shown in Figure 1 or explained in the text (e.g. Section 2.2). Add this info.

l34: « elasticity ». Delete word.

l34: « for retrieving sub-surface density and anisotropy ». In general, it is unclear to me if the parameterization is bijective, i.e., is there one unique anisotropy tensor and density for a given elasticity tensor?

l57: « Section 2 gives a theoretical overview of the elasticity tensor » -> « Section 2 gives the background of the elasticity theory ».

l.66: « Where the » -> « whose »

Eqn. 1: Give the assumption underlying this equation (Hill’s lemma).

l72-73: Explain what is « transversely isotropic » and that z is vertical (?).

Eqn. 2: Report also sigma and epsilon (as in Eqn. 1), so that the Voigt notation is explicit (there may be some variations with some 1/2, 2 coeff.).

l77: « common relations ». It would be convenient to have these relations in the appendix. Indeed, the paper change from one notation to another (C_ij, Lamé, bulk modulus, etc.) and it is sometimes difficult to follow.

l80-83: Only one Thomsen parameter is used after. Only present this one and explain in a few words what it represents. 

l94 : « 33 component » -> « the component C_33 »

l122: « elasiticity » -> « elasticity ». Check the orthograph in the whole paper with dedicated software to avoid typos.

l132: « HS bounds predict the effective elastic properties ». No, they are bounds (with one equal to zero).

l160: « influence of anisotropy increases monotonically ». Clarify if its relative anisotropy.

Fig. 1: show in log scale to be consistent with the rest of the paper. Show the expected range of validity of the models. « Illustration » -> « Evolution »

Table 1 : « Isothtermal » -> « Isothermal »

Section 3.4. Give reference to the choice of the ice properties.

Fig. 2: Are the first row and last column really necessary? You could gain space to make the subfigures larger.

Fig. 3: comparing C_FEM to C_G_33 (power law) is somehow unfair (scatter due to the fact that, e.g. C12 != C11). Indeed refits of the power law on each component show very little scatter (Tab. 2).

Fig. 6b: I am not sure this figure makes any sense. Anisotropy at high density affects elasticity anisotropy, but it appears that porous ice is not anisotropic (due to ice physics). See main comments. Can you make the same figure but with the FEM as the ground truth?

Sect. 4.6: This is not clear to me why epsilon_cryst should decay with increasing porosity. For sure, it cannot go above the value for a single mono-crystal. Moreover, you do not need this decay to draw your conclusion (geometric anisotropy is dominant for most of the densities). Simplify.

l287 : significantly

l340-342: You discuss here possible improvements. Does it really make sense with the given current performance and the uncertainty on the measurements? Delete paragraph?

l388: « The new parameterization constitutes a significant simplification ». I would not say it is simple but rather, « it is a crucial tool »

Reference :

• Hagenmuller, P., Matzl, M., Chambon, G., Schneebeli, M., 2016. Sensitivity of snow density and specific surface area measured by microtomography to different image processing algorithms. The Cryosphere 10, 1039–1054. https://doi.org/10.5194/tc-10-1039-2016
• Proksch, M., Löwe, H., Schneebeli, M., 2015. Density, specific surface area and correlation length of snow measured by high-resolution penetrometry. Journal of Geophysical Research: Earth Surface 120, 346–362. https://doi.org/10.1002/2014JF003266

Pascal Hagenmuller

(Please find attached a pdf version of the review. In the following maths symbols according to latex notations)

In this paper, the authors propose to use finite element simulations conducted on X-ray tomography images (395 images in total) to compute the homogenized elastic behavior of snow, firn and bubbly ice. The resulting behavior is modeled as transversely isotropic, which corresponds to 5 independent material parameters. Homogenizing the elastic properties of snow from X-ray tomography images is not new and several authors (cited in the paper) have already proposed such a procedure in the last decade. And some of them have already used a transversely isotropic model for snow. The contribution from the authors to the state of the art is to propose a fit over the whole range of porosity with the combination of a power law and the theoretical Hashin-Shtrikman bound (equations (9), (11) and (12)) to respect the fact the ice properties are recovered for a solid fraction of 1. The fit explicitly accounts for both density and geometrical anisotropy (estimated as the ratio of autocorrelation lengths in the vertical and horizontal directions). They show that their fit enable to achieve a higher precision than previous fit proposed in the literature.

Then, the authors discuss the relative contribution of geometrical anisotropy for different porosity values. The authors also assess the relative contribution of geometrical and crystallographic anisotropy on the elastic properties of snow, firn and bubbly ice. They show that the influence of anisotropy decreases with the decrease in porosity. They also show that geometrical anisotropy is dominant over crystallographic anisotropy up to a volume fraction of 0.7.

Even if the contribution to the state of the art is a little bit incremental on some aspects, I would suggest publication, provided the authors clarifies the following points. On the form, the paper is globally well written but the main story line is sometimes a little bit difficult to follow.

1. In Fig. 6. the authors explicitly show the relative contribution of geometrical anisotropy for different porosity. The authors could comment a little bit more this central Figure in their paper. For instance, there seems to be a tendency for $\alpha$ to increase with $\phi$ for low porosity values on the data set considered in Fig. 6.(a). Is there any physical explanation for that? In Fig. 6.(b) the two squares show that the larger over and under estimation zones are indeed not observed in the data set. Could the authors therefore comment on the maximum over and under estimations that one could get by not accounted for $\alpha$ for different snow densities? How does such uncertainties compare with uncertainties related to density estimations?
2. Time series of snow metamorphism are considered in the data base. In these time series (especially temperature gradient experiments), anisotropy develops. It could be interesting to show on some specific time series, how the fit propose by the authors enable to accurately capture the anisotropic evolution of the mechanical properties.
3. Section 4.3 may possibly benefit from some clarifications. I understand that Kohnen parametrization is valid at high ice volume fraction only. This could be stated explicitly in section 2.2.2. Then, why not having presented the results in the same form as in Fig. 3 with correlations between the different models and the FEM predictions?
4. I understand that the anisotropy is accounted in the $\boldsymbol{P}^\mathrm{ice}$ tensor in equation (9) which is related to the Eshelby tensor $\boldsymbol{S}$ recalled in Appendix A that depends the ratio $\alpha$ between the vertical and horizontal corelation lengths. Therefore, I do not understand why the tensors $\boldsymbol{M}$ and $\boldsymbol{M}*$ are introduced in section 3.3...
5. More details on the FEM simulations should be given. For instance, what are the boundary conditions?
6. In Fig. 2, when confronting the predictions of FEM against the U model, it could be nice to display the 1:1 line as done in Fig. 3. For the right graphs, the units (GPa) should be corrected as dimensionless quantities are plotted. Can the authors give more explicitly what is the expression of the fit curve? Does it refer to one of the specific models presented before? Interpreting the data in terms of Young or Bulk moduli could ease the physical interpretation of the parametrization. Instead, the authors simply refer to Torquato (2002a) to find the equivalences with respect to the coefficients $C_{ij}$.
7. In table 2, the formal expressions for the different models could be recalled or at least the number of the corresponding equations in the paper.
8. Fig. 5 is not very clear and do not bring much added value compared with Fig. 3... From Fig. 3 the authors have proved that their fit perform better than the other models. Why not using this depth profile to highlight the impact of accounting for the anisotropy or not in the PW model?
9. The data from Wautier et al. (2015) where snow is modeled with the same transverse isotropic behavior is available in the supporting information. Correlation lengths are also given. Maybe the authors could consider testing their fit on these data points? 

In the paper under review, the authors characterize the contribution of geometric anisotropy on elastic moduli for snow, firn, and bubbly ice. Specifically, the authors propose a normalized upper Hashin-Shtrikman bound for elastic moduli that encompasses a range in porosity from 0 to 1. Under this scheme, the geometric anisotropy ratio and fabric tensor is related to the elastic moduli by using an Eshelby tensor. The behavior of the elastic moduli is simulated using finite element methods via volume averaging on 395 images taken with X-ray tomography. Although the methods presented herein are not novel (as indicated in the referenced models), the normalization scheme presented in the present work provides an excellent fit to the simulated outputs for elasticity of dilute dispersion of spherical cavities and is, relatively, computationally inexpensive. Moreover, all five components of a transversely isotropic elastic modulus for snow, firn, and bubbly ice used in the present work were predicted using 2 parameters rather than 5 (required for simulations referenced in the present paper) in calculating an orthortopic elasticity tensor. The authors note the influence of both geometric and crystallographic anisotropy in the range of densities from snow to ice. At lower porosities, the contribution of geometric anisotropy is greater than that of crystallographic with a volume fraction around 0.7 (and has appreciable contribution to the elasticity moduli even at densities past the bubble close off density for firn/ice). At higher porosities, the influence of these two effects switches, such that crystallographic anisotropy dominates the behavior of the elastic modulus at greater depths in the firn/ice. However, the point at which this transition occurs is not resolved in the present study and should be a discussion for future work. Although not entirely novel, the present study provides the cryospheric sciences a new method for characterizing the elastic moduli across the range of porosities for snow to ice and the relative contributions of geometric and crystallographic anisotropy across the full porosity range (0 to 1 for snow to ice, as defined in the present study). I would suggest publication, provided the authors resolve the following major and minor points.

Major Comments

1. The authors remark in Sec. 2.2.1 that empirical parameters in Eq.4. need to be estimated by fitting to experimental data. At least it should be explicitly stated that constraints on a_ij and b_ij have not been made, and there has yet to be a widely agreed upon model based on laboratory and/or field measurements of snow to ice porosities. This is a serious limitation in comparing model outputs in the present work for the elastic moduli to that of the FEM simulations (and other model comparisons, such as presented on in Fig. 3.). Moreover, it would be nice to get a brief description on the conditions under which the a_33 and b_33 components were obtained.

2. In Sec. 3.2, the two-point correlation function is defined and computed via fast Fourier transform of the 3D tomography images. It is noted that, if using the model presented in Eq.(9), (11), and (12), which showed the best agreement to FEM simulations of elastic moduli compared to other models presented on in the study, only two parameters are needed to determine all five components of the elasticity tensor, ζ and β. A possible limitation of using an anisotropy parameter as defined in Sec 3.2. and the accompanying appendix, is that it requires knowledge of the correlation lengths using 3D X-ray tomography images, which may not be widely available or accessible to those in the broader snow, firn, and ice communities.

3. The authors should consider including the temperature time series presented on in Figures 6 an 7 and discussed in the concluding remarks. To that end, it is not clear, at least from how the model in Sec. 2.6 and the accompanying Appendix are presented, how the elastic modulus (or, similarly, Eshelby tensor) depend on temperature. It is clear that there is an effect on anisotropy that is due to temperature effects, however, without a formulation for the dependence of the anisotropy ratio or Hashin-Shtrikman upper bound of the effective elastic modulus on temperature one would expect it difficult to implement the model presented on the in the current work.

4. It should be noted early on in the present study how you are defining porosity and the reference frame you are using.

5. In equation 8, it is assumed that the dependence on the eigenvalues for the ice volume fraction are of power-law type. Why? One can ad hoc assume the relation follows a power law, but a more detailed explanation should be provided.

6. In figure 3(b), all components of the elastic modulus from the FEM simulations are compared to the C33 components of the power law model presented in Eq. 4. It may be beneficial to clarify why the density power law agrees more with C44 components (rather than the C33 for which other comparisons are made) obtained from the FEM simulations.

7. It would be nice to see a plot of the upper HS bound with the polarization or fabric tensor (as Srivastava et al. (2016) notes, the choice in which one does not effect the representation of geometric anisotropy).

8. Please provide a more explicit relation for effective elastic moduli (presented as Eq. (2) in the original text) to Young’s modulus, bulk modulus, and Poisson’s ratio. To that end, it would be useful to to see these relations plotted as a function of mass fraction for all discussed models.

9. On line 148, it would be useful to see a figure of the geometric result of α > 1, α < 1, and α = 1, to illustrate the result of prolate inclusions, oblate inclusions, and isotropic bubble distributions. Better yet if a movie of this transition could be provided across a range of porosities.

Minor Comments

I3: "... geometrical) that give rise to macroscopically anisotropic elastic behavior." to "... geometrical), which can give rise to elastic behavior due to macroscopic anisotropy."

I16: "...the elastic modulus is the probably..." to "the elastic modulus can be used to represent the mechanical properties of snow, firn, or bubbly ice, and so knowledge of the effective elasticity tensor plays a crucial role in..."

I22: "In particular,...anisotropy..." to "In Schlegel et al. (2019), the role of elastic anisotropy was emphasized. Specifically, the retrieval of elasticity..."

I24-26: "...anisotropic, on one hand...orientation" to "an anisotropic with respect to ice matrix geometry (e.g. ...) and crystallographic orientation [there needs to be a citation here]."

I29: "Recent work wave propagation..." to "Recent work by Hellmann et al. (2021) on measuring wave propagation in glacier ice suggests that at low porosity [give value] the effective elastic... is influenced by geometric effects (such as porosity)." Reduce the intensives (e.g. "already"). They weaken your argument.

I37: "Using the Finite-Element (FE) methods..." to "Using Finite-Element (FE) methods via volume averaging, a solution for static linear elasticity yields the material effective elastic properties."

I48-49: "the HS bounds incorporate the non-linear interplay between structural anisotropy and density." HS bounds incorporate the non-linear relation between density and bulk and shear stress, but you need to be more careful defining how anisotropy is represented in the limit of these bounds in the introductory remarks (or refer to the description in Sec.2.4).

I51-56: This entire paragraph is one sentence. Although this is fine, consider breaking it up to make your points more clear to readers.

I62: "by comparing it with the above mentioned shortcomings of previous work" to "... with previous work in which these parameters are not captured," or something similar. Refrain from adding subjective words.

I69: "... is defined by Hook's law" to "is defined by Hook's law, using Hill's lemma,..." Add the reference frame.

Eq.1. Consider adding the region over which the continuum is occupied. Also, consider adding a remark on the use of the notation in eq.1. in connecting volume averaged strain energy of a heterogeneous material at micro length-scales to that of a macroscopically heterogeneous material under uniform strain.

I86: Specify how you are defining ice volume fraction (see Major comments for a related remark).

I86: "...parameterization often..." to "parameterizations use a power law..."

Eq.6. Consider showing the limit explicitly.

I128-129: It can be left to the reader to refer back to the cited text. However, to avoid ambiguity, consider providing a brief description on how these parameters were obtained.

I132: "Hashin-Shtrikman..." to "when using Hashin-Shtrikman (HS) bounds, the effective elastic properties of porous materials can be derived based on volume fraction and microstructural anisotropy (incorporated through n-point correlation functions)."

I130: Consider re-phrasing the subsection header to specify the case of geometric(?) anisotropy, since the distinction is clearly made on I140-141.

Eq.10. Consider expanding on eq. 10 with 8. Also, make reference to accompanying definitions given in the appendix.

I155: "...the formulations including anisotropy, three different anisotropy ratios..." to "including geometric anisotropy, three different anisotropy ratios (alpha = 0.1, 1, and 1.6) were evaluated..."

I161: "...tend to an isotropic state" to "the geometric fabric must tend to an isotropic state"

I161: "the U bound" to "the upper bound (C^U)."

In Fig.1. Gpa should be GPa

Eq.11. Define the normalized HS bound before introducing the transformation.

Eq.13. This assumes no mass exchange between the two phases, correct? If so, please note this. I.e. that you assume no sublimation (a process that occurs in glaciological contexts and is a deviation in model applicability to natural environments).

I213: Consider including the definition of Q(\alpha) used here, for completeness.

I216-217: What temperature is this valid for? It would be useful to run FE simulations over a range of shear and bulk moduli that correspond to a range of temperatures. (See major comments)

I221-222: Consider adding a table summarizing the model equations, names, unknown parameters, and porosity range over which they are valid, and their main assumptions.

In Fig 3. Refer back to table 1. It is unclear in Figures 2 or 3 what the legend means. These codes should also be explained in the body of the text in Sec. 4.1.

In Fig 3.b. it is clear this is the only value for the elastic modulus for which the empirical parameters (a_33 and b_33) are known. However, either consider omitting this plot (3.b), since it adds confusion as to which models provide the best agreement to simulations of C, or obtain empirical parameters for a_ij and b_ij from other experimental datasets. Also, please explain a possible reason that C_33 from the density power law agrees more with C_44 from the simulations.

I247: Possibly refer to eq.(5) here.

I248-249: "...and with the literature P-wave velocity of ice...". Please include the conditions under which this was measured.

In Fig 4. "Gpa" to "GPa"

I258: "...for reasons discussed in Fig. 1" to "..., as mentioned in the caption of figure 1,"

I259: "...gives the right prediction" to "...parameterization provides an elastic modulus that agrees well with simulated values, taken from images of EastGRIP samples that were close to the surface. At these depths, anisotropy values are low (\alpha < 1), and are consistent with..."

I261:"... for deeper snow" to "at greater depths, the geometric anisotropy increases."

I262:"...demonstrate a good performance..." I am not sure what you mean by this, or at the least it is slightly vague. Please explain what you mean by good performance here.

I265: "... reasons discussed in Fig. 1" Omit and consider stating clearly the reason for greater error values for the ZC model with greater densities. Also, reasons cannot be discussed in figures (it is up to you to discuss what the figures mean). Please re-phrase.

In Figure 5: "...Bottom: Error plot which is given by the difference between the simulated elastic modulus..." to "Bottom: Error in FEM and PAR parameterized elastic modulus calculated from the difference between simulated elastic modulus..."

I274: "...with zero relative error by" to "with zero relative error for isotropic..."

I278: "...we show the geometrical Thomsen parameter \epsilon_geom (see Eq.3) in Fig. 7." Referencing this figure does not assess the geometrical vs crystallographic anisotropy in your calculations, or at least it is not clear what you mean by this. Consider "... we plot the geometrical Thomsen parameter, obtained from Eq. 3, against the porosity, the output of which is given in Figure 7."

I295-296: No parentheses are needed around the authors names.

I299: "In contrast... are explicit formulas." to "In contrast, the limiting behavior of the Hashin-Shtrikman bounds can provide an explicit formula for effective moduli."

In Fig. 7. Make sure the symbols you are using are consistent. For example, varepsilon is used to label the vertical axis, but epsilon is used in the text.

I303: "...collapsing onto" to "...collapse onto"

I304: "...as a function of normalized HS upper..." to "... as a function of the normalized HS upper bound..."

I305: "...parameterization Srivastava et al..." to "parameterization used in..."

I317: "...yields a eigenvalue zero..." to "yields zero eigenvalue..."

I323: "...where MIL resulted circle with no signatures of anisotropy." It is not clear what you mean by this. Please explain the outputs from the cited text more clearly. For example, "...observed by Klatt et al. (2017), in which a Boolean model of MIL with arbitrary rank fabric tensors, produced figures of circles when evaluated on Reuleaux triangles. Moreover, the MIL analysis used in that model was insufficient in detecting interfacial anisotropy."

I327: "...overall our parameterization shows" to "overall, the parameterization used in the present work (C_ij^PW, given by Eq. (11)-(12)), had excellent agreement (R² = 0.99) when fit to all components..."

I335: "... evident for temperature gradient time series (TGM2 and MMTO17) from Fig. 6 (a)...". Which one is MMTO17? This timeseries is not listed in the referenced figure (or, at least, it is not clear which time series you mean). Also, without a description of the temperature timeseries in the body of this work, the dependence of the elastic modulus used in the PW model (and even \alpha) on temperature is not entirely clear, other than vertically oriented structures being favored at high temperature gradients.

I371-375: "in principle... However, typos..." These sentences are not needed. Unless you plan to also compute crystallographic fabric at low porosity (such as in future work), with remedied typos from the referenced text, it detracts from the overall discussion.

This study derives a new parameterization for the effective elasticity tensor that is valid for the full range of volume fractions (i.e. for snow, firn, and bubbly ice). The authors compare this new parameterization to existing parameterizations valid for certain ranges of volume fractions, and identify the potential importance of geometrical anisotropy (in comparison to density and crystallographic anisotropy) in controlling elasticity.

The science and methodology appears sound and the results are interesting. My main comments are about the presentation of the material; in some cases I found the takeaways and the specific novel contributions of the work difficult to pull out of the descriptions. I would also recommend more description of some specific methods (possibly at the cost of some of the background material, which is quite extensive). Besides these recommendations, I would support publication.

Organization and Presentation

In general, I found the balance between the background/“literature review” section of the paper and the methods/results to be a bit off – there was quite a bit of background information, which in some cases was useful (it is helpful to know where the individual models come from and what assumptions they include) but the length and amount of information made it difficult to parse what novel contribution this study was providing. Further, as discussed further below, the background seems to come at the cost of description of methods, which I believe are important to include.

A section in the paper or an appendix that discusses what it takes to apply or use this tensor would be helpful. Similarly, I was left with questions about how generalizable this tensor is – the authors do a good job of explaining its generality in terms of volume fraction, but because the parameterization is based on empirically-found parameters, I believe it would be helpful to know two things:

• What are the conditions that these parameters are found in? What sizes of samples, grain sizes, temperatures, etc.?
• How well will this tensor generalize to different temperatures, grain sizes, etc.?

This would be useful in knowing how to apply this new parameterization.

Methods

I believe the paper would benefit from more detailed outlines of the methods used, particularly with respect to the X-ray tomography (how are the samples found/made? What conditions are they made/found in?) and the FE simulations (what is the resolution of the simulations? What are the assumptions underlying these simulations). Similarly, it would be helpful to have more information about the EGRIP samples – what is the specific variable identified in these samples (crystallographic anisotropy?).

Other Comments

• It would be potentially helpful to clearly define “geometrical anisotropy” up front before using the term. It is an important concept for the paper and for some audiences (including myself, since I do not study porous materials) the term is not obvious
• Figure 1: it would be helpful to include a more descriptive legend to remind the readers which tensor is meant to be valid for which ranges of volume fraction
• Equation 12: what is beta?
• Table 1: it would be helpful to have another column that included the region that each sample was obtained from (or if it was a laboratory sample)
• Figure 2: Why are there two legends? What is the difference between a-I and j-r? This information would be helpful in the caption
• Line 255: For some reason, I struggled to parse the sentence “Another view…our data”, which seemed important to understand what Figure 5 is showing.