Mapping 3D Structure of Loose Quaternary Deposits Combining Deep Learning and Multiplepoint Statistics: An example in Chencun, Northern Pearl River Delta
 ^{1}School of Earth Sciences and Engineering, Sun YatSen University, Guangzhou, 510275, China
 ^{2}Guangdong Provincial Key Lab of Geodynamics and Geohazards, Guangzhou, 510275, China
 ^{3}Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai), Zhuhai, 519080, China
 ^{4}Guangdong Geology Survey Institute, Guangzhou, 510080, China
 ^{1}School of Earth Sciences and Engineering, Sun YatSen University, Guangzhou, 510275, China
 ^{2}Guangdong Provincial Key Lab of Geodynamics and Geohazards, Guangzhou, 510275, China
 ^{3}Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai), Zhuhai, 519080, China
 ^{4}Guangdong Geology Survey Institute, Guangzhou, 510080, China
Abstract. Reconstruction and cognition of structures of the Quaternary deposits, like thickness variation and displacement, is necessary for understanding neotectonics and the evolution of palaeovalleys and deltas. Multiplepoint statistics (MPS) is a useful method to reconstruct threedimensional geological models in many fields. However, nonstationary spatial patterns and semantics in geological blocks are difficult to extract and reconstruct with the MPSbased methods, especially for those probabilitybased MPS methods. To reconstruct 3D characteristics of loose Quaternary deposits and the semantic relationship between them, an algorithm coupled MPS and deep artificial neural network (DANN) is proposed. The DANN is constructed and used to extract and simulate the global characteristics of geological structures. Process of sequential simulation and stratigraphic sequence calibration are implemented to build an initial model. To obtain a reasonable final realization, an iterative MPS simulation process with a multiscale strategy is implemented. With several crosssections and trench profiles used as modeling dataset, two concrete examples of constructing the Quaternary sediments in Chencun, South China are given. The displacements of sedimentary formation belonging to the Pleistocene reveal the strata rupture caused by the fault activities. The modeling results illustrated that the DANN used in the method can extract and simulate global structures of Quaternary deposits, and MPS simulation with the ExpectationMaximizationlike iteration process can optimize local characteristics in results effectively.
 Preprint
(3321 KB) 
Supplement
(0.42 KB)  BibTeX
 EndNote
Weisheng Hou et al.
Status: final response (author comments only)

RC1: 'Comment on gmd202283', Anonymous Referee #1, 21 Jul 2022
This paper gives a new solution for the model construction with DANN and MPS. The DANN extract global characteristics of geological structures. and the MPS simulation with ExpectationMaximizationlike iteration process for local characteristics of the geobodies.The idea is interesting and may attract the focus of the modeller.Because in real reservoir,to guarantee the long range conenctivity and characteristics of small architectural elements is often a difficult task for modellers. The method will be adopted for modeling facies in oilfield.In my opinion the paper is excellent and should be accepted if some comments are addressed:
1.The cross section for global structure is parallel,why? is it sufficient for the reproduction of the characteristics .In another direction,the structure may be different, how to reproduce in 3D space with DANN? please clarify
2.The small structue is reproduced with the constrain of trenches to update the initial global structure. is it a 3D template for comparison.How to guarantee the consistency between the global stuctures and the local trenches. May be an illustration is needed.
3. In line 261263. the method of constructing 3D models from 2D training image isstudiedï¼So maybe a revised depiction is suitable.

AC1: 'Reply on RC1', Henggaung Liu, 08 Aug 2022
Thanks for the comments.
 The cross section for global structure is parallel,why? is it sufficient for the reproduction of the characteristics .In another direction,the structure may be different, how to reproduce in 3D space with DANN? please clarify
Response:
In many cases, 3D model is constructed with parallel crosssections. Therefore, we provide a real case with parallel crosssections in this study. However, as you pointed out, the parallel crosssections cannot sufficient information to reproduce structures in the other direction. Therefore, as shown in the example in section 4.2, seven nonparallel trenches are used as data source for modeling. In the proposed method, whether the crosssections used for modeling are parallel or not is not a key issue.
In this study, the modeling area is partitioned into regular grids. Then, the core of modeling is to specify the geological attribute on each grid. To simplify the problem, geological objects are formed with closed surfaces, which the same attribute appear in it. So, the problem is transferred into obtaining the elevation of subsurface as the interpolation process.
Here, the objective function of the DANN is constructed with the elevation of subsurface. Note that each object has two networks for top surface Mmax（i） and bottom surface Mmin（i）. Values of the top surface and bottom surface are normalized and used as training dataset to learn the spatial distribution of geological surfaces. Here, coordinates (x, y) of the simulation grid are the input data, and the corresponding elevations hmax(x,y,Atti) and hmin(x,y,Atti) of each geological attribute are labeled for training Mmax（i） and Mmin（i）.
The DANN is trained with the depths of contact lines from the crosssections. Then, the elevations of the top and bottom surface of the Att_{i} are predicted when the coordinates of grids to be simulated are input into the Mmax（i）and Mmin（i） . After completing the traversal, the top and bottom surface of each attribute can be obtained.
 The small structue is reproduced with the constrain of trenches to update the initial global structure. is it a 3D template for comparison.How to guarantee the consistency between the global stuctures and the local trenches. May be an illustration is needed.
Response:
The two models are independently built by the presented method, although they are in the same area (shown in the followed figure). The model with trenches are not constrained by the global model. Two trenches of BTI2 and BTI3 crosses the TI3 for the major model. In these two cross positions, the geological contacts are consistent in the trenches and parallel crosssections. Because the TIs and BTIs are hard data for the modeling, the geometries of the two model in the cross positions are kept consistency.
 In line 261263. the method of constructing 3D models from 2D training image isstudiedï¼ŒSo maybe a revised depiction is suitable.
Response:
As pointed out in lines 261263, the 2D TIs cannot provide information in third dimension. In this study, to obtain 3D patterns for modeling, we provided a compromise way by expanding 2D images to 3D data rather than directly extract 2D patterns for constructing 3D structures. The detail process of expansion is provided in lines 265286. To express the process clearly, here, a simple process is given as follow:
 The first step is defining an expansion area Buff_{sec }with or bigger than a template size in the simulation grid (SG), in which the 2D TI is in the middle or on the boundary. The position of Buff_{sec}depends on the location of 2D TI.
 A grid node neighboring to the known data or to the node with appointed data is randomly selected and marked as the current access grid nodeu_{c}_{ }in the lateral layer. The attribute that appears with the maximum number in the moving window of 3×3 grid nodes of which the center node is u_{c} is assigned to the attribute value of u_{c}.
 Repeating step (3) until to values of all grids are assigned inBuff_{sec }layer by layer, the 3D TD with a template size is obtained.
Note that the presented method is implemented with multiple scales. Therefore, in each scale, the 3D TD should be constructed before simulation.

AC1: 'Reply on RC1', Henggaung Liu, 08 Aug 2022

RC2: 'Comment on gmd202283', Anonymous Referee #2, 08 Aug 2022
The paper combines two methodologies, MPS and DL, to create realistic geological models in a very specific case study.
A number of claims are made regarding literature that are not correct, plus the literature review is significantly lacking in major recent contributions in this area
 “the MPS method difficult to reconstruct global spatial features with anisotropic and nonstationary characteristics”. Many methods of nonstationary MPS exist, please consult the chapter “Nonstationary MPS” in the book of Mariethoz & Caers, 2015.
 “However, threedimensional structure cannot be directly extracted from 2D crosssections in the MPSbased simulation method”. Many methods exist that use 2D cross sections to create: e.g. Comunian, A., Renard, P. and Straubhaar, J., 2012. 3D multiplepoint statistics simulation using 2D training images. Computers & Geosciences, 40, pp.4965. This problem (of stereology) is quite common in 3D imaging (e.g. Xray; MRI) and many methods of DL exists to do this. CNNs are very popular for this.
There is also no mention of the recent contribution of GAN methods starting with Laloy
Laloy, E., Hérault, R., Jacques, D. and Linde, N., 2018. Trainingâimage based geostatistical inversion using a spatial generative adversarial neural network. Water Resources Research, 54(1), pp.381406.
Song, S., Mukerji, T., Hou, J., Zhang, D. and Lyu, X., 2022. GANSimâ3D for conditional geomodelling: theory and field application. Water Resources Research, p.e2021WR031865.
In addition, one would wonder if the case study specified would not be better solved with surfacebased or implicit domain geological modeling. This literature is also missing:
Frank, T., Tertois, A.L. and Mallet, J.L., 2007. 3Dreconstruction of complex geological interfaces from irregularly distributed and noisy point data. Computers & Geosciences, 33(7), pp.932943.
Consider an example with much larger complexity than presented in this manuscript:
Yang, L., AchtzigerZupanÄiÄ, P. and Caers, J., 2021. 3D modeling of largescale geological structures by linear combinations of implicit functions: Application to a large banded iron formation. Natural Resources Research, 30(5), pp.31393163.
Also consider rulebased geological modeling methods: Pyrcz, M.J., Sech, R.P., Covault, J.A., Willis, B.J., Sylvester, Z., Sun, T. and Garner, D., 2015. Stratigraphic rulebased reservoir modeling. Bulletin of Canadian Petroleum Geology, 63(4), pp.287303.
The question for me is: what is the new contribution of this paper? What has been achieved that is new and hence exportable to other cases, applications? My take on this question is that the contribution remains narrow
 The methodology seems tailored to their specific case study. As a result, it contains many, many adhoc choices and tuning parameters that I would not know how they extend to other cases.
 The methodology is directly applied to the case study, there is no other verification, for example, we do not learn about how it would apply to some simpler synthetic models.
 Because of the many adhoc tuning, the methodology is very complex for what may be easier solved with other geological modeling approaches. There is small likelihood that others would use this method for that reason. The paper does not contain any comparisons, except for broad methological comparison which (see above) are not always accurate.
The manuscript is a significant amount of work, and a lot of thinking went into modeling this specific case study. But then, my question for the editor is if this is sufficient for publication in GMD where the aim would be to share modeling approaches across many applications.

AC2: 'Reply on RC2', Henggaung Liu, 15 Sep 2022
Dear Editor,
We appreciated the referees’ valuable comments on the manuscript titled with “Mapping 3D Structure of Loose Quaternary Deposits Combining Deep Learning and Multiplepoint Statistics: An example in Chencun, Northern Pearl River Delta” (gmd202283). We revised the manuscript in accordance with the comments. We hope the manuscript will meet your standards for the next process.
All revised text is marked in red in the manuscript.
Sincerely,
Weisheng Hou and Hengguang Liu et al.,
School of Earth Sciences and Engineering
Sun YatSen University, China
The following are onebyone responses to the referee’s comments. The comments are shown in black. Our replies are in red, with parts of the text copied from the revised manuscript for your convenience.
The paper combines two methodologies, MPS and DL, to create realistic geological models in a very specific case study.
A number of claims are made regarding literature that are not correct, plus the literature review is significantly lacking in major recent contributions in this area
Response:
We appreciated that the referee pointed out the defects in the manuscript. Also, the text is revised carefully and some references are added to the revised manuscript.
 “the MPS method difficult to reconstruct global spatial features with anisotropic and nonstationary characteristics”. Many methods of nonstationary MPS exist, please consult the chapter “Nonstationary MPS” in the book of Mariethoz & Caers, 2015.
Response:
The statement here is not clear. We do agree that many MPS approaches have been proposed to simulate 3D geological structures with nonstationary characteristics as Mariethoz & Caers did (2015). Here, we addressed breaking through the limitation of conditional probability calculation and the local optimization in the MPS stochastic simulation with constraints of global features. Because in the process of MPS simulation, the continuity of geological blocks is dependent on the candidate patterns and their relationships. In the scenario of constructing layered strata or more complex geological structures like a sedimentary cycle, intrusive rock, and overthrust structures, constraints with global constraints become vital. The existing methods for simulating nonstationary structures, such as auxiliary variables or images, are not fit for the scenario mentioned above. For example, the Herten aquifer structure is constructed with the MPS method with constraints of the Kriging interpolation (Coumian et al., 2012). From this, we gave this statement in the text. We deleted this statement in the revised text to clear up the misunderstanding.
 “However, threedimensional structure cannot be directly extracted from 2D crosssections in the MPSbased simulation method”. Many methods exist that use 2D cross sections to create: e.g. Comunian, A., Renard, P. and Straubhaar, J., 2012. 3D multiplepoint statistics simulation using 2D training images. Computers & Geosciences, 40, pp.4965. This problem (of stereology) is quite common in 3D imaging (e.g. Xray; MRI) and many methods of DL exists to do this. CNNs are very popular for this.
Response:
Here, the sentence is not exactly stated. Many MPS methods such DS (Mariethoz et al., 2010), s2Dcd (Comunian, A., Renard, P., Straubhaar, J., 2012. 3D multiplepoint statistics simulation using 2D training images. Computers & Geosciences 40, 4965.), CCSIM algorithm (Tahmasebi, P., Hezarkhani, A., Sahimi, M., 2012. Multiplepoint geostatistical modeling based on the crosscorrelation functions. Computational Geosciences 16, 779797.), CCSIMTSS method (Ji, L., Lin, M., Jiang, W., Wu, C., 2017. An Improved Method for Reconstructing the Digital Core Model of Heterogeneous Porous Media. Transport in Porous Media 121, 389406.) and SNESIM (Strebelle, S., 2002. Conditional Simulation of Complex Geological Structures Using MultiplePoint Statistics. Mathematical Geology 34, 121.) can build 3D models by using 2D crosssections directly. Some of these methods extract the conditional probability distribution from single or multiple twodimensional TIs or generate 3D models by staking a series of twodimensional slices. Also, the hybrid method (Gueting, N., Caers, J., Comunian, A., Vanderborght, J., Englert, A., 2018. Reconstruction of ThreeDimensional Aquifer Heterogeneity from TwoDimensional Geophysical Data. Mathematical Geosciences, 50, 53–75) generates multiple twodimensional slices as conditional data and TIs. We have discussed it in detail (Hou et al., 2021; Hou et al., 2022). We also observed that these methods did not use the 3D patterns directly in the simulation when 2D crosssections are used as TIs. So, getting 3D patterns from 2D crosssections directly is a good alternative option in patternbased MPS methods. Using the 3D patterns in calculating similarities between candidate patterns is more intuitive and can reduce the artifacts in stochastic simulation. We presented an interesting solution to construct 3D training data from 2D crosssections (Hou, W., Liu, H., Zheng, T., Shen, W. and Xiao, F., 2021. Hierarchical MPSbased threeDimensional geological structure reconstruction with twodimensional image(s). Journal of Earth Science, 32(2): 455467; Hou, W., Liu, H., Zheng, T., Chang, H., & Xiao, F. (2022). Extended GOSIM: MPSdriven simulation of 3D geological structure using 2D crosssections. Earth and Space Science, 9, e2021EA001801).
Therefore, the sentence is revised as:
Directly extracting 3D patterns from 2D crosssections should be concerned in the MPSbased simulation method (Hou et al., 2021, 2022). In this study, the 2D geological crosssections are extended into 3D TD. The expansion process is described as follows:
There is also no mention of the recent contribution of GAN methods starting with Laloy
Laloy, E., Hérault, R., Jacques, D. and Linde, N., 2018. Trainingimage based geostatistical inversion using a spatial generative adversarial neural network. Water Resources Research, 54(1), pp.381406.
Song, S., Mukerji, T., Hou, J., Zhang, D. and Lyu, X., 2022. GANSim3D for conditional geomodelling: theory and field application. Water Resources Research, p.e2021WR031865.
Response:
Thank you for your kind reminder. Some more DL methods have succeeded to build 3D geological models in the past several years. In the revised manuscript, we added some new contributions including these two papers. In addition, we also stated the difference between the artificial neural network we used and these algorithms in the revised manuscript. The data we use for training are known data with the incomplete research area, rather than existing 3D models that can be used as modeling to provide a 3D space model. Because in the case mentioned in this paper, it is difficult to obtain the 3D model that contains the study area's geological characteristics. Moreover, in this case, there are not enough 3D models to support the training of GANrelated generation models.
In addition, one would wonder if the case study specified would not be better solved with surfacebased or implicit domain geological modeling. This literature is also missing:
Frank, T., Tertois, A.L. and Mallet, J.L., 2007. 3Dreconstruction of complex geological interfaces from irregularly distributed and noisy point data. Computers & Geosciences, 33(7), pp.932943.
Consider an example with much larger complexity than presented in this manuscript:
Yang, L., AchtzigerZupanÄiÄ, P. and Caers, J., 2021. 3D modeling of largescale geological structures by linear combinations of implicit functions: Application to a large banded iron formation. Natural Resources Research, 30(5), pp.31393163.
Also consider rulebased geological modeling methods: Pyrcz, M.J., Sech, R.P., Covault, J.A., Willis, B.J., Sylvester, Z., Sun, T. and Garner, D., 2015. Stratigraphic rulebased reservoir modeling. Bulletin of Canadian Petroleum Geology, 63(4), pp.287303.
Response:
Thank you for your kind reminder. Indeed, the cases in this study can be realized by many other methods. Also, much large complexity could be built by rulebased methods or implicit methods.
Thank you for your kind reminder. Indeed, the cases in this study can be realized by many other methods. Also, much large complexity could be built by rulebased methods or implicit methods.
In this study, we presented a different idea to construct layered strata, faults, and intrusive rocks, which is difficult for the MPS simulation method.In the case used in this paper, the global spatial characteristics of these geological objects are shown as strong anisotropy, singularity of spatial distribution, continuity of geometric form of geological objects, such as the fracture zone. Moreover, compared with other articles mentioned, the algorithm used in this article's example is based on a small number of geologicial profiles, without a large number of borehole and cross sections as constraints. The data we used are sparse and can be uneven distribution. The study is more concerned with the method to get global features for modeling, and the way to introduce semantic constraints into the model process. At this stage, the proposed method really cannot handle scenarios with complex geological phenomena like reverse strata.
The question for me is: what is the new contribution of this paper? What has been achieved that is new and hence exportable to other cases, applications? My take on this question is that the contribution remains narrow
 The methodology seems tailored to their specific case study. As a result, it contains many, many adhoc choices and tuning parameters that I would not know how they extend to other cases.
 The methodology is directly applied to the case study, there is no other verification, for example, we do not learn about how it would apply to some simpler synthetic models.
 Because of the many adhoc tuning, the methodology is very complex for what may be easier solved with other geological modeling approaches. There is small likelihood that others would use this method for that reason. The paper does not contain any comparisons, except for broad methological comparison which (see above) are not always accurate.
The manuscript is a significant amount of work, and a lot of thinking went into modeling this specific case study. But then, my question for the editor is if this is sufficient for publication in GMD where the aim would be to share modeling approaches across many applications.
Response:
According to the loss function of the MPS simulation (Eq. 1), we cannot solve the function directly, because we only have TI.
(1)
Obviously, it is a typical underdetermined problem. An EMlike algorithm is a good alternative to solve this problem. Note that the EMlike algorithm needs an initial model for the following iterative process. Here, if an initial model with unreasonable structures is used, the final realization may not be accepted.
Therefore, in this study, the core idea is to create a 3D initial model with global features extracted from known data by deep ANN, and then an MPS simulation process is implemented on the initial model, where the semantics between geological blocks are introduced into the simulation process. We think the main contributions of the algorithm proposed in this paper are as follows：
 The deep ANN constructed in this paper can obtain the global spatial characteristics of geological objects with limited and unevenly distributed data, instead of using existing threedimensional models. The GAN or other DL methods usually need rich data for training the network. Based on the initial model constructed by the obtained global spatial features, the local spatial features of the simulated geological objects are optimized by using the multipoint statistics algorithm.It makes the final model more consistent with prior knowledge.
 The geological semantic is extracted from the crosssections directly and is used as a constraint in reconstructing sedimentary deposits. The algorithm obtains reasonable stratigraphic sequences from the training image, and these stratigraphic sequences are used to identify and adjust the attributes where stratigraphic sequence errors occurred in the threedimensional reconstruction with MPS simulation, which makes the 3D model in line with prior knowledge. In addition, the algorithm also combines the multiscale iterative process to reduce or even eliminate the artificial artifacts caused by the correction of stratigraphic sequence.
The algorithm constructed in this paper is universal, and the parameter comparison in this paper is only to determine the response characteristics of various parameters used in the algorithm to the simulation results. The algorithm comparison in this paper is mainly reflected in the comparison with the ExtendedGOSIM algorithm.
In addition, we have actually used different cases to verify the algorithm. Due to the length limitation of the article, we did not show these cases in the article. This is the reconstruction result of the algorithm applied to an example with completely different geological conditions (Fig. 1).
Figure 1. An engineering example application based on the algorithm proposed in this paper. (a)~(f) shows the distribution of each geological attribute in the simulation results. (g) shows the 6 geological profile data used for modeling.

CEC1: 'Comment on gmd202283', Juan Antonio Añel, 23 Aug 2022
Dear authors,
After checking your manuscript, it has come to our attention that it does not comply with our "Code and Data Policy".
https://www.geoscientificmodeldevelopment.net/policies/code_and_data_policy.html
You have archived your code on Mega. However, Mega is not a suitable repository according to our policy. Therefore, please, publish your code in one of the appropriate repositories (e.g. Zenodo, PANGAEA, FigShare,...) and reply to this comment with the relevant information (link and DOI) as soon as possible, as it should be available for the Discussions stage. In fact, with the problems listed in this comment, the manuscript should have never been accepted for Discussions.
Also, you use a methodology based on Deep Learning and Artificial Neural Networks. Therefore, you must include in the repository the primary input and output data used in your work. In the same vein, we do not accept that you share binary files, and a binary file does not let to check and evaluate the code and precludes reproducibility. Therefore, please, you must upload the source code and any instructions needed to reproduce your experiment and run your model.
Also, I have not seen a license listed in the Mega repository. If you do not include a license, the code continues to be your property and can not be used by others, despite any statement on being free to use. Therefore, when uploading the model's code to the repository, you could want to choose a free software/opensource (FLOSS) license. We recommend the GPLv3. You only need to include the file 'https://www.gnu.org/licenses/gpl3.0.txt' as LICENSE.txt with your code. Also, you can choose other options: GPLv2, Apache License, MIT License, etc.
In this way, you must include the modified 'Code and Data Availability' section in a potential reviewed version of your manuscript, the DOI of the code (and another DOI for the dataset if necessary).
Be aware that failing to comply with this request will result in the rejection of your manuscript.
Juan A. Añel
Geosci. Model Dev. Executive EditorJuan A. AñelGeosci. Model Dev. Executive Editor
AC3: 'Reply on CEC1', Henggaung Liu, 15 Sep 2022
Dear Editor,
We appreciated the referees’ valuable comments on the manuscript titled with “Mapping 3D Structure of Loose Quaternary Deposits Combining Deep Learning and Multiplepoint Statistics: An example in Chencun, Northern Pearl River Delta” (gmd202283). Thank you for your suggestions and reminders. We have replied to all reviewers and will modify the article, and upload the code and case we built again as required. But it will take some time. Thank you again for your advice and patience.
Sincerely,
Weisheng Hou and Hengguang Liu et al.,
School of Earth Sciences and Engineering
Sun YatSen University, China

AC3: 'Reply on CEC1', Henggaung Liu, 15 Sep 2022
Weisheng Hou et al.
Weisheng Hou et al.
Viewed
HTML  XML  Total  Supplement  BibTeX  EndNote  

402  90  19  511  21  4  2 
 HTML: 402
 PDF: 90
 XML: 19
 Total: 511
 Supplement: 21
 BibTeX: 4
 EndNote: 2
Viewed (geographical distribution)
Country  #  Views  % 

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