Response to Referee #2 (Florian Wellmann)

> In the manuscript with the title “Blockworlds 0.1.0: A demonstration of

> anti-aliased geophysics for probabilistic inversions of implicit and

> kinematic geological models”, the authors present an approach to overcome

> the aliasing effect of forward-calculated geophysical fields from non-ideal

> subsurface mesh discretizations. The approach and the provided implementation

> in a geological and geophysical modeling framework is relevant, as it enables

> an efficient treatment of such aliasing effects, and also timely, as these

> types of subsurface parameterisations are widely used in geophysical inverse

> schemes, especially in cases where regular and rectilinear meshes are used.

> The approach will therefore be of use in many similar modeling cases.

> The manuscript is very well written and follows the logical structure of a

> scientific document. The introduction is quite detailed, but provides a

> comprehensive overview over the current approaches at the link between

> geological modeling and geophysical inversions and, in my point of view,

> this detail is also important to place the developments into the frame of

> current developments. The authors provide detailed descriptions and

> references to similar approaches. The methods section provides a very

> instructive motivating example, before delving into the detail of the

> implementation. Even if this order may seem a bit odd, it is common practice

> in computer science publications and here a good choice to quickly show what

> the problem is and how it can be mitigated, before describing the approach

> in detail.

We thank the referee for their comments, which align well with our intentions in writing the paper.

> The antialiasing method is presented with kinematic modeling examples,

> where the interaction of multiple geological events quickly leads to complex

> structural settings. Even if this is suitable to evaluate the use of the

> method, it is not directly evident that it will also work in other cases

> where the block model is obtained through different interpolation approaches,

> for example. The authors explain the current state of modeling very well in

> the introduction - but especially implicit methods, mentioned even in the

> title, are not evaluated in detail. I would suggest to clarify this point

> and also to remove implicit models from the title.

We respectfully disagree with Referee #2 in regard to the evaluation of our technique with implicit models.  Our view is that while Blockworlds is a kinematic model, it is also an implicit model, in which the positions of interfaces are determined by values of a scalar field.  The mathematical implementation carries over directly to other interpolation methods like co-kriging, and thus to models such as GemPy or LoopStructural, as long as the scalar field and its gradient can be evaluated within each mesh cell.  In fact both co-kriging and kinematic models exist for many real regions so it isn't a question of representational power.

We can't see any reasons why our method isn't general enough to carry over to other interpolation methods that involve level sets of scalar fields, and so we will defend the use of "implicit" in our paper's title.  However, if Referee #2's concern is more about the range of geological structures we've tested (as with Referee #1), we have now included new model runs using fold configurations that target specific hypothesized shortcomings in our method.  We have also added some text to the Introduction and the Discussion to clarify and circumscribe the relation of Blockworlds to other kinds of models, moving the more general formulation covering other interpolations to Section 5.1 (new version lines 472-491).

> More minor aspects that could also be clarified or commented on:

> 

> 1.  What is the influence of data error? Would the aliasing effect be

>     similarly relevant or would it decrease?

Increasing the data error makes the likelihood less restrictive, which does relax the spurious fragmentation of the posterior surface by aliasing effects.  This can be viewed as a form of tempering.  However, it does not change the terracing effect along the directions of spatially aliased parameters mentioned in comment #2 below.  We see similar effects when decreasing the sample size.

In our spherical intrusion example we find that for Gaussian data errors, we have to inflate the noise variance by a factor of 5-10 -- so that the data error per point is of the order of the full dynamic range of the gravity signal -- before the different modes merge.  We have included a new figure in section 2.2 (now Figure 2) that shows this.

> 2.  The problem of block representations is also that small changes in

>     the parameters lead to no change in the likelihood function, and

>     therefore to problems in specific inverse algorithms (not only in MCMC).

>     It would be interesting to explain this aspect in more detail.

We emphasize MCMC in our paper because we see the capacity for full posterior inference as one of the major advantages of parametric geological inversions.  Optimization algorithms seeking a single point estimate, for example maximum likelihood or maximum a posteriori, could easily miss the effects we describe.

It is true nevertheless that the parameter bias caused by aliasing would still affect even point optimization algorithms.  Gradient-free methods, such as downhill simplex, may seem to work if data constraints are weak but are liable to get stuck in local modes.  Ensemble methods such as differential evolution might do better at finding multiple aliased modes.  Gradient-based methods will in general fail badly because the likelihood is piecewise constant along the directions of aliased geometric parameters, so trajectories of dynamics using an aliased gradient will be poorly matched to the larger-scale posterior geometry.

We have added a new paragraph along these lines to Section 2.2 (new version lines 163-175) which introduces the aliasing problem.  If there are particular other methods that Referee #2 would like us to comment on, we would be happy to do so.

> 3.  There are several geophysical inversion schemes that are not based on

>     “blocky” mesh structure, but on pillar grids that conform better to the

>     exact interface location (e.g. IGMAS+

>     https://igmas.git-pages.gfz-potsdam.de/igmas-pages/)- but are then

>     limited in terms of there possibility to represent geologically complex

>     models. But it should be mentioned that the (vertical) aliasing-effect

>     would be very limited for these methods. 

We appreciate the pointer to IGMAS, which seems to use unstructured tetrahedral meshes for solving the finite-volume potential field problem.  This has elements in common with the Cauchy surface approach of Zhdanov and Liu (2013) that we cite in our Conclusion, which indeed should experience no aliasing effects as we have described them here.

Referring to Wellmann & Caumon (2018) for a description of "pillar grids", we agree that we expect a large reduction in aliasing effects in the vertical direction, though the benefit would be lower for interfaces or faults cutting across the pillar axis.

We have added this discussion and a citation to IGMAS to our Conclusions (new version lines 530-535).

> 4.  The anti-alaising approach works when the exact location of the

>     interface inside a cell can be determined. Many geophysical inversion

>     schemes first perform the mapping to a fixed grid, and then perform the

>     inversion (I believe this is the case for the cited approaches by

>     Giraud,etc.). In this case, the anti-aliasing approach would not work,

>     correct? It would be good to mention this aspect (maybe also in the

>     limitations section in the discussion).

The objective of the Giraud et al. inversions is different, as we point out in our Introduction; it uses an assumed geological block model as the prior mean in a voxelised inversion for rock properties.  No inference is done directly over the geological parameters, as in Wellmann et al. (2018), de la Varga & Wellmann (2018), Olierook et al. (2021) and the present work.

In the Giraud et al. case, if the projected rock properties are aliased, then the prior mean is biased.  It is still possible to adjust the rock properties away from the prior mean to improve the fit to the data, and the likelihood in the voxelized representation should remain continuous.  Anti-aliasing the block model used as the prior mean should improve results relative to an aliased prior, though whether the improvement is significant may depend on the model error term (see our response to Referee #1's comment about model error).

Distinction between the Giraud workflow and our workflow is essential to the interpretation of what we are doing.  We draw this distinction in the Introduction, but we now reiterate it in Methods section 2.2 (new version lines 183-191).

Response to Referee #1 (anonymous)

> The manuscript presents a method for addressing the problem of defining a

> geological model on a regular grid and its effects on forward as well as

> inverse geophysical modeling. The problem being that structural geological

> models consisting of planes and surfaces loose definition when exported to

> a regular grid for geophysical modeling, and consequentially that this has

> an effect on physical modeling. This problem is termed an “aliasing problem”.

> The new contribution in this work lies the construction of an anti-aliasing

> function that captures the forward geophysical response more accurately than

> a regular grid discretization. The anti-aliasing function is then applied in

> the context of structural models construction using kinematic models. Using

> a synthetic case study, the effect of aliasing, and hence anti-aliasing is

> studied.

> Discretization of continuous geological models, in particular structural

> models is an important problem and the anti-aliasing method is an

> interesting idea, that has merit for a publication.

We thank the referee for their comments, and are pleased to hear that they see merit in our approach to undersampled interfaces in discretized geology.

> The manuscript suffers however from unfocused writing and presentation.

> The main contribution is the anti-aliasing function. There is no new

> methodology in kinematic modeling and Bayesian inverse modeling using McMC.

> It is unclear to me why the Authors decided to focus specifically on

> kinematic models, or on a rather complex example of McMC. It would be

> sufficient to study the effect of anti-aliasing on the misfit function

> or likelihood model. To illustrate their proposals, a series of models

> (not just from kinematics) could have been constructed, and the method

> tested, without the need to perform inversions. Instead, examples and

> effect of the anti-aliasing function are described mathematically,

> but not using examples.

We see our contribution differently.  The likelihood arises in the context of an inference problem, describing the probability of the data given the model parameters.  We're not sure it makes sense to specify a likelihood outside the context of an inference problem, that is to say, an inversion.  Without varying model parameters, the most we could do would be to examine the distribution of residuals for fixed geology, and we agree that would be a very different paper from the one we have written.

The way the inference is set up also matters.  If the rock properties are themselves the model parameters, the likelihood will still be a continuous function of its parameters even if an undersampled block model is used as the reference model (prior mean).  The rock properties could be directly adjusted to fit the data, unless the model error term was very restrictive.  On the other hand, in parametric inversions that condition rock properties directly on structural parameters -- and thus, implicitly, use a very restrictive model error term -- you see the effects we are describing.  See our discussion of model error below, and a new paragraph in Section 2.2 (new version lines 183-191) that draws this distinction explicitly.

We set all this up in the Introduction, which, though quite detailed, is meant to situate our work within the context of the current developments linking geological modeling with geophysical inversions.  We think the detail is justified given the broad readership of GMD.  Although the elements of what we are doing are not new in themselves, we have not come across a thorough description of the problem we solve nor a synthesis of these elements in the recent literature we have read (and cite) in this area.  We developed the anti-aliasing prescription because we saw this as the minimal modification

to whole classes of geological modeling tools necessary to enable them to be used widely within *parametric* inversion workflows.  And we are interested in parametric inversion workflows because these enable full posterior inference (using MCMC) in ways not easily done for flexible voxelized models.  Our Introduction thus seems, to us, not much longer than needed to explain the impact we see for our work.

We focus on kinematic models because we found kinematic models to be easy to implement from scratch as a demonstration.  We found existing model codes to require much more up-front effort to patch for people who were not the authors (and harder to convince the authors to implement without some demonstration of a new feature's value).  But many real areas have both kinematic models and other 3-D models (say, from co-kriging), so the classes of structures Blockworlds can represent are comparably rich.  We would not consider adaptive Metropolis to be a very complex MCMC algorithm; rather, we think it is the simplest algorithm that gives reasonable performance without excessive hand-tuning on the badly scaled posteriors that arise in these problems.

Throughout specific edits mentioned below we do our best to clarify two key points:  that anti-aliasing solves a problem found in integrated geology-geophysics workflows that may not seem obvious when examining the elements of such workflows on their own; and that the rock properties do not, in our current formulation, appear as freely adjustable parameters, but are deterministically mapped from the geological parameters.

> I suggest therefore that the manuscript be rewritten with the anti-aliasing

> method be more thoroughly studied and evaluated on examples that better

> show the strength and limitation of the anti-aliasing function, outside 

> the context of inversion first, then only apply inversion as an illustration

> case, but in my view this is not needed. This would also mean considerably

> shortening the introduction, which discusses all kind of matter tangential

> to the anti-aliasing problem and have more to do with inverse modeling,

> which is not the main contribution of this work.

We understand the value of examples, and so have added several new models intended to stress-test the elements we highlight as potential limitations; limitations are now discussed in Section 2.4 (new version lines 234-254).  Our first batch of models covered faults, so this second batch covers folds (new version lines 304-320), in light of the "banded iron formation" flagged in subsequent comments.

We have also reworked the Introduction, shortening it somewhat and reordering material to improve the flow.  However, it contains substantially the same range of citations and content, which as explained we find necessary to explain what we see as the problem and why we chose to solve it in this way.

> The presentation of the main idea, section 2.3 is unfocussed, it takes

> 3 paragraphs of meandering thoughts about inverse modeling to get to the

> main point/contribution. The presentation then is not very clear; it is

> never explained what the reasoning is behind Eq (2)-(4). Then suddenly,

> line 205, a “training appears”, yet there is not yet any motivation as to

> why this is done. My suggestion is to write more structured and intentional,

> first explaining the intent of what is going to be done and why, then only

> the how. From what I understand, the anti-aliasing function is a proxy

> regression models that is fitted on a dataset created within a single voxel.

> It would be better to say this first.

We have reworked the first half of Section 2.3 in line with this feedback (new version lines 201-214).  We find the text to be vastly improved and thank Referee #1 for the helpful suggestions.

> ** Please refrain discussing from what “might” be done, unless evidence is

>    shown. Line 218

We now omit this line.  In light of this and subsequent comments, we have also gone through the manuscript and removed other instances of vague or conditional language, clearly labeling elements as future work where applicable.

> ** Line 220: “Should”, is it or isn’t’ it? Otherwise please omit.

Omitted.

> ** Paragraph starting at 221: this whole section reads as non-committing,

>    do you in fact have solution for voxels of varying size or for

>    curvi-linear coordinates, or is this mere conjecture?

Now clearly labeled as future work.

> ** Anti-aliasing for implicit models. This looks interesting, but I believe

>    it is no longer studied later, so this would not be a contribution unless

>    it is fully worked out and illustrated.

We have moved this text (new version lines 472-491) to the Discussion, which we have renamed "Discussion and Future Work".  We emphasize in that section however that the mathematics carries over to any model that uses isosurfaces of scalar fields to define interface positions, of which Blockworlds is one.

> ** Line 243: how is Eq (6) an “expansion”. If it is an expansion, what are

>    the higher order terms?

Changed to "first-order approximation".

> ** Where does Eq (9) come from?

We now explain this in more detail.  The idea is that the anti-aliasing function provides a smooth transition in rock properties at the interface.  The way the function was developed in Section 2.3 ensures that it gives the correct answer when evaluated at the grid centers.

> ** Line 250: “it should be accessible to implicit models based on

>    co-kriging, and could be spot-checked only for cells lying near

>    interfaces”. It should or it is? This very terse sentence needs

>    expanding, what is the role of co-kriging here? If your sentences

>    contain “should” and “could” then perhaps they should be omitted,

>    or moved to discussion/future work.

On further reflection in responding to this point, we realize that the construction to which this is referring (the Hessian quadratic form) is not a very useful metric.  Without a calibrated correction to the anti-aliasing function, this just says anti-aliasing works better for low-level curvature, which we already knew.  It also can't reliably flag actual undersampling of geological structures in the Nyquist sense if it is evaluated on the same mesh as those undersampled structures.

We have therefore cut material on these metrics and simply mention these limitations at the beginning of Section 2.4 (new version lines 234-240).

Instead, we suggest the Kullback-Liebler divergence as a single number summarizing the correspondence between versions of the posterior calculated with different mesh resolutions and anti-aliasing settings (new version lines 396-402).  Unlike the Hessian-based measures, we actually implement this and add the values as an additional diagnostic to the MCMC metrics summary (now Table 3).  We stress that this is a global metric on probability distributions that is most easily calculated by MCMC for these problems, supporting the use of anti-aliasing in the setting for which we developed it.

> ** Line 253, what maximization is carried out? What is the function being

>    maximized?

Now omitted (see above).

> ** Where does (10) come from?

If Referee #1 means Eq. (10), this is earlier in the text.  We agree that the definitions of terms in these expressions are incomplete, and have revised Section 3.1 (new version lines 271-286) to try to make this clearer.

> ** Line 256: true but a sphere has the simplest possible curvature

> ** Line 261: a folded surface has no curvature?

Both omitted.

> ** The section 243 – 265 discusses the possibly main issues with the

>    contribution, in my view needs to be worked out with real examples

>    of folding with different variation in curvature and how the authors

>    plan to address it. Can you anti-alias a complex banded iron formation?

>    Banded iron formation (BIF) - Stock Image - C029/7116 - Science Photo

>    Library

We believe Referee #1 is referring to Section 5.3, which covers lines 433-469 in our first version.  The question of where the curvature breaks down is important and is reasonable to address now instead of in future work. We have moved this entire section to Section 2.4 (new version lines 234-254) to reflect its importance.

As mentioned above, we have added additional models with folds of different wavelengths, some with folded folds or folded faults, to explore the behavior of layers with curvature and complex structure.

If the description of a banded iron formation comes from an implicit model, then we certainly can anti-alias it.  However, this configuration seems irregularly folded in ways that are challenging to model outside of co-kriging or a similar framework.  Since Blockworlds is meant to be an accessible demonstration at this stage rather than a full-featured modeling package, we feel it is reasonable to proceed with a parametric investigation of adversarial cases as suggested above.

> ** Line 258 – 273: this needs extensive examples and illustrations

> ** Section 3.1 needs figures to explain parameters and relations to geometry

We intended for Figure 3 to provide the examples and illustrations.  If Referee #1 means for us to label the diagrams with parameter names, that should be relatively straightforward.  We have accordingly shortened Figure 3 to show examples of kinematic parameters for different events, and have replaced Figure 4 to show the final versions of various models.

 

> I will not be reviewing the next section on inverse modeling in detail,

> since these are not new contributions. My recommendation is to work out

> better the anti-aliasing idea, illustrated with more cases, understanding

> the extent of applicablity, then resubmitted, possibly omitting inversion.

See discussion above.  If Referee #1 has a specific kinematic model in mind that uses an event history composed of event types Blockworlds supports, we will be happy to run it.

> A final comment lies in understanding the amount of approximation relative

> to other errors. To me what is presented seems relevant when you have very

> well constrained problems, e.g. gravity inversion with simple geological

> structures. Once you increase the complexity of the models, then the

> posterior uncertainty increases as well, hence the approximation become

> much less relevant. In addition, one would also need to compare the error

> discretizing with data and model error. Why would you need to anti-alias

> when there are uncertainties orders of magnitude larger? Knowing this

> would be a considerable contribution.

Referee #2 also raised the point of data error.  We have taken a quick look at the amount of data error needed to merge the modes in the sphere demonstration in Section 2.2 (new version lines 176-182), and find that we need to inflate the data error by a factor of 5-10 -- making errors much larger than residuals -- or to reduce the number of data points from 400 to

16.  Increasing data error or deleting data points also does not change the terracing effect in which the derivatives of the likelihood with respect to aliased parameters is zero except on a discrete set of discontinuities.  It thus does not address the potential failure of gradient-based algorithms.

Our current statistical model assumes zero model error, resulting in a direct projection of the geological parameters into the data space.  This is not common practice in geophysics, but as near as we can tell is consistent with what has been done in other papers with similar workflows, such as Beardsmore et al. 2016 and Wellmann et al. 2018.  We feel an important future direction would be to use parametric prior means in voxelized inversions, as mentioned in our Conclusions -- producing models that can both fit the data flexibly and perform geological inference by marginalizing over the rock property variables.  In this scenario treatment of model error is also important.

For nonzero Gaussian data and model error, the voxelized rock properties can be integrated out in closed form, resulting in a modified Gaussian data likelihood with strong correlations.  It is therefore easy to imagine that sufficiently large model errors will temper or relax the posterior in ways similar to larger data errors, but the behavior as a function of model error

will be more complicated than for data error.  We tried to implement this, but were unable to make a numerically stable version to demonstrate by the response deadline.  We agree that the issue is important, but beyond the scope of this paper -- given the potential computational challenges, and since we have shown that aliasing effects arise even for otherwise perfectly specified models with weak data constraints and fine voxelizations.  We hope to revisit model error in a later follow-up paper.

The usual model error term in voxelized geophysical inversions addresses one facet of the broader question of (parametric) model mis-specification, which is important but also beyond the scope of our paper to treat comprehensively.

We believe aliasing will actually become a bigger problem as the complexity of the model increases, especially in kinematic models, because the rich conditional dependence structures caused by overprinting will create vastly more complex posterior geometries.  We hope to address this in subsequent papers using the machinery we have built here.

Response to Chief Editor (Juan Antonio Anel)

> After checking your manuscript, it has come to our attention that it

> does not comply with our Code and Data Policy.

> https://www.geoscientific-model-development.net/policies/code_and_data_policy.html

> You have archived your code in GitHub. However, GitHub is not a suitable

> repository. GitHub itself instructs authors to use other alternatives for

> long-term archival and publishing, such as Zenodo. Therefore, please,

> publish your code in one of the appropriate repositories, and include

> the relevant primary input/output data. In this way, you must include in

> a potential reviewed version of your manuscript the modified 'Code and

> Data Availability' section, the DOI of the code (and another DOI for the

> dataset if necessary).

We thank the Executive Editor for bringing this to our attention.  The Topical Editor also mentioned this before the start of the review process, at which time we acquired a Zenodo DOI (10.5281/zenodo.5195426, dated Aug 13) reflecting the state of the code when producing the plots and tables as seen by the Referees.  We communicated this DOI but forgot to change the GitHub URL in the manuscript itself.

We have updated our code in response to the Referees' feedback, and have included the updated release (DOI: 10.5281/zenodo.5759225, dated Dec 5) in the Code and Data Availability section.  Since the datasets we use are synthetic data generated deterministically by the code covered under this DOI, we don't believe we need a separate DOI to cover explicit copies of them.