Underworld-GT Applied to Guangdong, a Tool to ... - Springer Link

Journal of Earth Science, Vol. 26, No. 1, p. 078–088, February 2015 Printed in China DOI: 10.1007/s12583-015-0517-z

ISSN 1674-487X

Underworld-GT Applied to Guangdong, a Tool to Explore the Geothermal Potential of the Crust Steve Quenette*1, Yufei Xi2, John Mansour1, Louis Moresi3, David Abramson4 1. Monash eResearch Centre, Monash University, Clayton 3800, Australia 2. The Institute of Hydrogeology and Environmental Geology, Chinese Academy of Geological Sciences, Shijiazhuang 050061, China 3. School of Earth Sciences, The University of Melbourne, Melbourne 3010, Australia 4. Research Computing Centre, The University of Queensland, St Lucia 4072, Australia ABSTRACT: Geothermal energy potential is usually discussed in the context of conventional or engineered systems and at the scale of an individual reservoir. Whereas exploration for conventional reservoirs has been relatively easy, with expressions of resource found close to or even at the surface, exploration for non-conventional systems relies on temperature inherently increasing with depth and searching for favourable geological environments that maximise this increase. To utilitise the information we do have, we often assimilate available exploration data with models that capture the physics of the dominant underlying processes. Here, we discuss computational modelling approaches to exploration at a regional or crust scale, with application to geothermal reservoirs within basins or systems of basins. Target reservoirs have (at least) appropriate temperature, permeability and are at accessible depths. We discuss the software development approach that leads to effective use of the tool Underworld. We explore its role in the process of modelling, understanding computational error, importing and exporting geological knowledge as applied to the geological system underpinning the Guangdong Province, China. KEY WORDS: Underworld-GT, geothermal potential, computational modelling approach, Guangdong. 0 INTRODUCTION Knowledge of heat flow and temperature in the shallow crust is the basis for estimating a reservoir’s geothermal energy potential. A component of resource estimation is often to model heat flow anomalies within a target region, searching for high temperature at shallow depth as a primary control on potential. Uncertainty of observational data is another dominant factor in estimating resource potential. Local heat-flux/depth observations generally have well-understood error characteristics and hence may be leveraged to better constrain models. However, as the exploration approaches depths beyond those readily accessible to direct observation of the thermal gradient, it is useful to consider constraints and interpretations at the crust and mantle scale that have been obtained by analysis of multiple datasets to be proxies to direct observation. Subsequently the complexity of the physical model and the data used to describe heat flow may evolve through experimentation whilst exploring. We can contrast this with historical approaches to heat flow oriented exploration, for example, extrapolated temperatures and heat flow measurements (e.g., Somerville et al., 1994) and Kriging-based 1D models (e.g., Chopra and Holgate, 2005). These methods maximise the utility of the direct observational *Corresponding author: [email protected] © China University of Geosciences and Springer-Verlag Berlin Heidelberg 2015 Manuscript received June 20, 2014. Manuscript accepted October 8, 2014.

information (mostly from shallow drill-holes)—surface heat flow measurements and bulk heat production and conductivity data, but no way to incorporate other available knowledge of the region, for example depth and lateral heterogeneity in the geology. The source data tend to be sparse and of variable quality yielding under-constrained results (Musson et al., 2009). Moreover, these methods and their tools are difficult to extend to incorporate any other constraints. Whether the driver is resource exploration or fundamental questions in geoscience, agility in the descriptions of heat flow and of the constraints it is of paramount importance. The aim of Underworld-GT (Quenette and Moresi, 2010) is to leverage the benefits of the Underworld tool (Moresi et al., 2007) into regional scale geothermal problems. Underworld is designed to permit the following to readily change throughout experimentation: (a) the models that describe the physics (primarily, the heat flow equations) and (b) the models that describe how data is used to constrain results. Described here is one of a four part collection of works that utilise novel geophysical and software methods to describe the heat flow of Guangdong. Descriptions of the inverse methods employed to achieve formally constrained geological/heat-flow models will be described in another paper. Similarly the application of the techniques to explain heat flow in the Guangdong Province will be described in another paper. 1 APPROACH The following is an abstraction of how models described by software, change as a result of the modelling process, con-

Quenette, S., Xi, Y. F., Mansour, J., et al., 2015. Underworld-GT Applied to Guangdong, a Tool to Explore the Geothermal Potential of the Crust. Journal of Earth Science, 26(1): 78–88. doi:10.1007/s12583-015-0517-z

Underworld-GT Applied to Guangdong, a Tool to Explore the Geothermal Potential of the Crust sidering both evolving physics and evolving data assimilation (Quenette, 2014). Let us consider the forward model form : →

(1)

where represents a set-theoretic map of elements in the domain of reals to the elements in the co-domain of reals . In effect, denotes the number of inputs to and denotes the number of outputs, implying and as vectors when the cardinality is greater than 1. Here we will assume only real numbers ( ) for now, but clearly, in software, values can take many other forms. In this form, the function may represent, for example, the heat flow equation (see Equation 12 below) implemented over some domain, or it may extend to mean a whole tool such as underworld implementing heat flow (with potentially many models working in unison, such as finite element discretisation, and so on). For simple heat flow, the cardinal could be the number of nodes in the domain where is the temperature at each node. The cardinal could be an enumeration of material properties (such as conductivity κ) for each distinct unit of material in the domain. If we consider there is an element of and for each location in the geology, as we delve deeper it is increasingly likely parts of and/or are not known. Employing a finite element discretisation, for example, to solve for across the domain will produce a relatively large . The input parameters over the domain tend to be abstracted into some relatively smaller number ( ′ ) of distinct material units. We generally expect ≫ ′ for to always be true. To understand heat flow at depth the intention is to use forward models (an instance of ) as the predictor over geological space where we otherwise do not have information. Hence we are intrinsically interested in the inverse function of the form : → (2) Once again we expect significantly lower numbers ( ) of geological observables ( , the input or known parts of , for example known temperature throughout the domain), and hence ≫ . Subsequently ≪ and thus is non-unique. There are many ways to address such inverse problems that utilise the forward model , rather than , to approximate the unknown parts of and . A potential property of this approach is that ′ . If we assume a Bayesian approach we might also view in the form ′ : → m+o=n (3) ,

′

′

,

where denotes the unknown parts of , denotes the known parts (parameters) of , denotes the known parts (observations) of , ′ is the unknown parts (predicted) of , and is the observation noise. Note that the co-domain of ′ is over both the known and unknown parts of . As such the fitness of any realisation is given by the error between as observables and as predicted by . Changes in the combination of data used in and , for

79

example the addition of new observations, implies a different number of models enacted in unison to form ′ , for example additional terms may be added to the model equation. The knowledge of which data/models is required to meet the goal are not known in advance. Hence part of the forward model itself can be considered to be unknown, such that we may view ′ as a function of the form : ℙ → Domain ∈ :

→

Domain π ∈ ℙ , where ℙ ⊆ ′

,

(4)

π, ,

Here we introduce algebra by the domain as a generalisation of functions ( ). Functions take some numbers and produces some other numbers. A sub-model to our total problem takes the form of domain ℙ, a subset of , meaning such sub models can also be represented as a function. The total intended forward model ′ is now some base-line model that additionally takes some further models π as input. The relationship that needs to hold then is ′

≡

∘π

(5)

That is, the decomposed version ( ∘ π) is equivalent to the total function intended ′ ). If we consider common programming languages, such as C and Python, they are implementations of , providing the user a language to define new functions ( s). A user can for example, define a function that takes some input real numbers and produces a real. However, these languages don’t readily implement the composition operator (∘) described in (5). Such a programming language keyword would take two functions and create the effect of a combined one, whilst ensuring the equivalence. Yet if such an operator existed, the equations and algorithms that define and solve a model problem could be written as software modules cleanly encapsulating those features. It would be easier to explore new equations. Further, existing optimisers would work just as well as if the software was written for performance. Some further consideration and annotations are required to ∘ compositions to ensure rigour in the equivalence. A key concept is to constrain ∘ to be a feature compositor. In software engineering, the concept of variability is defined as the tendency of something in the real world (which may or may not be a concern of the software) to change that is brought about by intent (that is not by accident) (Pohl et al., 2005). Object oriented programming is an enabler of variability within source code. However the objects themselves only exist at run-time, and consequently their use eliminate many important compiler optimisations. Feature oriented programming is an approach where each model requirement (a feature) has a direct one-toone mapping from design to source. In programmer terms— each model becomes its own ‘object’ and the use of this object in the code is programmed outside of any other model object. However, these ‘objects’ exist at compile-time. In AHEAD (Batory et al., 2003) notation, an application (i.e., code) can be described by a hierarchical collection of feature relationships, for example ,

,

(6)

Steve Quenette, Yufei Xi, John Mansour, Louis Moresi and David Abramson

80 ,

ψ

(7)

∘

∘

,

∘

,

(8)

where and are an encapsulation of feature requirement and respectively. They can have some arbitrary number of parts available to variation. In this example , and could be behaviour, some state structure and some other behaviour respectively. The feature is the composition (∘) of and parts, for some compositional operation relevant to each part. The approach in modelling problems is to treat π as chosen variants within a family of scientific codes. For example let us assume the feature requirement is the constitutive behaviour (9) where is heat flux, is temperature, and is thermal conductivity of the material. For now, (9) forms part of a chosen steady-state heat flow model, the feature , given by ∙

(10)

where is volumetric heat production. Hence the total feature represented by (10), which we call , is given by ∘

(11)

Which yields effect of ∙

(12)

But where (9) and (10) are considered and coded as two isolated feature concerns, bound some time later by the AHEAD-oriented instruction given in (11). It is considered that the feature crosscuts (Kiczales et al., 1997) both models—it has a mixture of both the heat flow model and the constitutive behaviour. Typical implementations of lexically compose these two features (that is, written as one function written in source code because its simpler to read and encourages compile-time optimisations). But it does not accommodate variability well—that is when (9) or (10) change. Furthering the example, let us assume is the Finite Element method, and that the behaviour is an element matrix assembly routine. The finite element assembly of the heat flow model ( ∘ ∘ ∘ ) is some appropriate composition of those three features. In this approach and are programmed within the encapsulations of , and so on (as separate feature objects). That is to code only one assembly routine and maintain each model that constitutes the stiffness matrix as lexically orthogonal. Further, the binding of a given model to the assembly is also described orthogonally to the assembly, at some phase after the model code has been written but before the code is run. A programming approach that implements composition, as an orthogonal sequencing of execution, is aspect oriented programming (AOP) (Wand et al., 2004; Kiczales et al., 1997). StGermain (Quenette et al., 2007) is a greatly simplified implementation of AOP for feature-oriented composition (where the qualities of expressibility and obliviousness are relaxed, and the static property is lacking). It approximates the aforementioned approach for dynamic scientific codes. For example, the StGermain dictionary of components is an approximation of the

AHEAD approach. Continuing the example, let us extend the heat flow to consider temperature dependent conductivity. In Danis et al. (2012) the following example is used (13) for . Here and denote the thermal conductivity at surface temperature and at some critical temperature respectively, beyond which thermal conductivity is held constant. is the surface temperature. The implementation of this non-linear constitutive behaviour is as another feature layer as a cross cutting concern, such that ∘ ∘ , where is the equation (13) impact on the finite element assembly. However, equation (13) impacts other concerns, such as the material model that now needs , , and to accommodate the inclusion of this non-linear solution approach. 2 UNDERWORLD The tool Underworld-GT is a collection of feature increments on Underworld, which in turn is a collection of increments upon StGermain. The geothermal toolbox within Underworld-GT includes some tools for prescribing inverse geological heat-flow problems and analysing them. Underworld is a 3D, parallel, geodynamics-modelling software toolbox. It provides governing equations, discretisations, constitutive relationships and analysis required for geodynamics model building and experimentation. It has been designed to scale from laptops to 1 000s of distributed-memory parallel computation (Quenette and Hodkinson, 2009). Underworld builds upon MPI (Gropp and Lusk, 1996) and PETSc (Balay et al., 1997). Underworld also includes gLucifer (Stegman et al., 2008), a toolbox to Underworld that computes visualisations in tandem with the solver, thus producing visualisations “on the fly” leveraging the finite element interpolation operations where ever possible. Some applications of Underworld include: coupled models of subduction with an upper plate (Fig. 1), mechanical evolution of the lithosphere (Fig. 2) and heat flow in systems of basins (Fig. 3). A Eulerian finite element discretisation underlies most

Figure 1. Subduction of an oceanic plate under a continent leads to the formation of an Andean orogeny and a concave trench morphology. Similar dynamic coupling is invoked to explain the topography of the Andes as well as the Bolivian Orocline, above the subducting oceanic Nazca plate. In grey: contour of the viscosity field, topography in color is interpreted from the vertical stress, arrows: velocity field (Capitanio et al., 2011).

Underworld-GT Applied to Guangdong, a Tool to Explore the Geothermal Potential of the Crust

81

We use Radial Basis Function (RBF) interpolation to map the material properties at the sample points of the structural model onto properties of the integration points. Leveraging our approach and software, the features of RBF and PIC are encapsulated and composed in a systematic and reusable manner. The RBF techniques employed are similar to those applied in the Smoothed Particle Hydrodynamics (Monaghan, 1992). For example, let the function define a geological surface at arbitrary location in the domain. , a smoothed interpretation of is given by ′

Figure 2. Shear bands forming during simple-shear deformation of a layered system of frictional material overlying a viscous substrate.

′

′ where , finite support, and

d

′

,

d ′

(14)

is a truncated Gaussian-like kernel with

1

(15)

The parameter determines the support size, and therefore also the degree of smoothing. We form a discrete approximation to equation 14 ̅≡∑

Figure 3. Modelled heat flow of the Sydney, Gunnedah and Bowen basins of the east coast of Australia (Danis et al., 2012). In this model we propagate the lithological properties from the calibrated Sydney-Gunnedah thermal model as an initial model into a green field. The 150 ºC isotherm is plotted in cyan. Topography is over laid. Major units shown. Underworld applications. Further numerical models augment this, either chosen in advanced or as part of the problem parameterisation. The Particle In Cell (PIC) method—a modified material point method (Sulsky et al., 1995) implemented in the Underworld toolbox PICellerator, substitutes the standard (Gaussian) quadrature module which is commonly used in most FE codes with a quadrature rule based on material point positions (Moresi et al., 2003). For the steady-state thermal model in which there is no deformation (12) we choose to keep the standard Gaussian quadrature. The properties of PIC are reused here as the basis for consistent, concise and robust discretisation of the heterogeneous geological material coefficients, whilst interpolating to where the finite element unknowns exist on a regular quadrilateral/hexahedral mesh. 3 IMPORTING We wish to populate the domain with material properties from (for example) geological interpretations of geophysical data, where the distribution of properties constitutes an assumption or describes one member in an ensemble realisation. This requires an accurate, efficient and robust mechanism for mapping arbitrary input data to these material points, which in-turn provides a consistent approach to data importation.

,

Δ

(16)

This equation allows us to determine an interpolation of the discrete data points, A ≡ , located at arbitrary locations in the domain. The infinitesimal volume element d ′ is replaced with the finite volume element Δ . This representation lifts the discrete data into continuous functions, which are subsequently applied to mean geological surfaces. Where surfaces are provided as point-sets, ̅ is similar to a signed level-set approximation of the surface. Where surfaces are provided as voxel-sets, nearest neighbour values can be used, or alternatively ̅ can be formed to provide a weighted average of neighbouring cells. These two methods cover geological models from tools such as GoCAD (Mallet et al., 1992, http://www.pdgm.com/products/gocad), Surfer (Lafaurie et al., 1994), and GeoModeller (http://www.geomodeller.com). 4 EXPORTING The StGermain plugin to Underworld (provided in the Underworld-GT toolbox)—DrillCores, implements virtual drill hole samples. The plugin itself is a component singleton that hosts an arbitrary number of drill holes. Each drill hole is described by a latitude and longitude (or for example northings and eastings), and some arbitrary number of equally spaced sample points down the hole. Either the depth of the hole (against a mean sea level (MSL) datum such as mAHD) and elevation (in m), or the total length (as measured depth (MD)/ ground level (GL) datum) must be given. Note that using the ground level datum can be undesirable as the evaluation of the topography surface (ground level) varies with resolution. An inclination and azimuth can also be given. The drill core model then provides for each drill hole sample point , ,

,

,

, ,

,

,

, ,

,

(17) (18)

That is: the evaluation of temperature, conductivity and three components of heat flux via the finite element basis functions at each point. is a given drill hole and is a given


82

0.0

Conductivity (W/(m·K)) 1.0 1.5 2.0 2.5 3.0

0.5

4.0

-100

50

0

-50 (a)

-100 290

280

300 310 320 330 Temperature (K)

340

-50

100

Temperature Observed Conductivity

Depth (mAHD)

Depth (mAHD)

100

3.5

Heat flow (mW/m 2) 0 50

Heat flow (eastings) Heat flow (northings) Heat flow (vertical)

50

0

-50 (b)

-100 350

100

-1.0

-0.5 0.0 0.5 Heat flow (mW/m 2)

1.0

Figure 4. (a) The observed (Danis, 2014) vs a simulated geotherm for the GDB1 bore hole in the Sydney Basin. A transition in conductivity in the 100 m from 50 mAHD is also shown, producing an inflection in the simulated geotherm. (b) The lateral (top axis scale) and vertical components (bottom axis scale) of heat flow with depth. Note the complexity in heat flow signal, from otherwise low frequency changes in material properties in the geo-logical/physical model. observation point on that drill hole. The depth is in MSL. The heat flux is given by (18) or as defined by the end-user as operations on the available fields. An example of this is plotted in Fig. 4. The drill hole evaluation model shares a cross cutting concern with the finite element evaluation and the various fields in the domain. From an aspect-oriented perspective, we state that the drill hole evaluation is advised onto the solver, keeping the solver code free from concerns of drill holes and specific fields. 5 COMPUTATIONAL SCALING The discretisation, importing and exporting models enable a robust mechanism to scale computational fidelity of an experiment. The following weak-scaling model can be used to assess the computational error as an experiment parameter available for optimisation. ,

∙

,

∙ ∙∑

,

2∙ ∑

,

∙ , ,

,

∙

, , ,

(19) (20)

That is to weak-scale the computational domain as per (19), where is the number of mesh elements in each component direction ( , , ), is the scaling starting from 1, is the number of drill holes, is the number of observation points for a given drill hole and is the total number of observation points across all the drill holes. The error at a given scale is then the difference of each observation point of each drill hole at the current and prior scale. The error might also be evaluated using the bounds of min and max, rather than the mean. Hence we become primary focussed with increased resolution near observations. To this end, Underworld-GT provides a mesh refinement feature that ‘node-bunches’ to a particular region of interest. 6 GUANGDONG The geological province of Guangdong belongs to South China fold belt, and is located in the southeast of the Eurasian plate. The province is bound to the north by the Yangtze Craton

Figure 5. A model of the Cooper Basin (Meixner and Holgate, 2009), where here we employ the mesh refinement to concentrate the element density in the region of the drill hole data. The mesh deforms to the topography also. and to the west by the Youjiang Block (see Fig. 7). It boasts the third largest number of active hot springs in China, including 7 springs over 90 centigrade (see Fig. 8 for the hot springs and boreholes considered here). These hot springs exhibit properties similar to conventional geothermal systems when observed to shallow depths. Yet this region is at some distance to plate margins or other obvious heat sources to conventional systems. The locations of hot springs also spans an unusually large lateral region, and the geological understanding of structure and heat flow at depth is relatively poor. Thin lithosphere and crust under Guangdong Province is observed in seismic tomography (Zhang et al., 2013; Zhou et al., 2012). Significant insight is provided through considering the tectonic setting. The present tectonic framework of Guangdong is controlled by the complex interactions between the Eurasian plate, Pacific Ocean plate, the Philippine Sea plate and the Indian plate. The Yanshanian movement (a counterclockwise rotation of Chinese continent, with westwards subduction and


Figure 6. Using models (18) and (19) on the Cooper Basin model to ascertain model sensitivity to resolution (scale). The horizontal axis is the total number of elements. The vertical axis is the error in degrees Kelvin. compression of Okhotsk and Izanagi plate (Wan, 2011)) during the Mesozoic formed the upper mantle and transition zone structure. Since the Late Mesozoic the Southeast China lithosphere and crust has thinned dramatically, resulting in the increase of deep thermal tectonic activity, asthenosphere upwel-

83

ling, and development of the coastal Cenozoic rift basin, volcanic eruptions and earthquakes (Ma, 1987). According to the reconstruction of (Honza and Fujioka, 2004), the Southeast China trench was located near the Guangdong coastline at the Late Cretaceous (~85 Ma) and gradually became inactive after the Early Oligocene (~35 Ma), transitioning to the Izu Bonin trench in the Late Oligocene (~25 Ma). Tomographic studies (Li and van der Hilst, 2010) indicate widespread low-velocity anomalies in the upper mantle at 100 to 400 km depth beneath the Guangdong Province. This anomaly may be an expression of wedge material from the prior Southeast China trench, likely significantly strained through the dynamic process to and the curvature of the nearby Izu Bonin trench. Mesozoic granitoid-volcanic rocks are distributed widely in Southeast China (Zhou et al., 2006). The instantaneous heat generated during magmatic activity affected the local geothermal field, and the radioactive heat producing elements in granite rocks provide both further increased temperature and delayed cooling of the crust in this region (Zhang et al., 2007). Further, our analysis of EGM 2008 and GRACE gravity data below suggests a step in the Moho parallel to the Guangdong coastline. However, further geometric heterogeneity of heat producing structures need to be explored, and can be achieved through a modeling process. Described here is a heat

Figure 7. The position of Guangdong. The background is topography (SRTM 30, Farr et al., 2007; USGS, 2004). Blue lines denote the main tectonic units: SB. Sichuan Basin, YB. Youjiang Block; black line: province boundaries; red line: Guangdong Province, study area; black line: South China and Japan trenches in the Late Cretaceous (~85 Ma) and Earliest Eocene (~52 Ma) as well as the Izu Bonin trench during the Latest Oligocene (~25 Ma), based on Li and van der Hilst (2010).

84

Steeve Quenette, Yufei Y Xi, John Mansour, M Louiss Moresi and Daavid Abramsonn

ngs and boreholes. LC. Lech hang; RH. Renhua; NX. Nanxxiong; SZ1. Sh henzhen1; SZ22. Figure 8. Thee position of naatural hot sprin Shenzhen2; WY. W Wengyuan n. a Grant, 19700) for the sedim ment layer, thee analysis (Spector and baseement and the Moho. M The topography is SRT TM 30 (Farr ett al., 2007; 2 USGS, 2004) 2 and surfaace temperaturee is the averagee clim mate temperaturee for the past 50 years (NCEP P). Whereas thee topo ography is givenn as a surface oof the model, accommodating a g a preesent limitationn in Underworldd the basal surfa face is flat at 355 km from f sea level and with the M Moho above it. The T model alsoo inclu udes a calculation of the Curiee temperature (5 578 ºC). Ta able 1 Examplee parameters tthat meet the Curie C surface constraint on the Guangd dong heat flow w model

Figure 9. A 3D D model map of o Southeast China C capturing the Guangdong Province, P facin ng the northeast. The top surrface is observed toopography and d sea level. Th he second surfaace is the top of thee basement. Th he third surfacce is the calcullated Curie depth. The lower moost surface is the Moho. Thee vertical axis has been exaggeraated by 10 timees. flow model off the Guangdong Province thatt provides a basseline for further invvestigations intto the geologiccal properties of o the province. We use a 3-layer model, which includees the major geological units of sediment, s basem ment and mantle (see Fig. 9).. The region spans 740 7 km×4 000 km×35 k km, axis aligned to zonne 38. The EGM 20008 model, including sphericaal harmonic cooefficients to the degree and orderr of 2 159 (Pavlis et al., 2008), was used in conjuunction with thhe GRACE (grravity recoveryy and climate experiiment) satellite mission. We seeparated gravityy into depths using multi-scale waavelet decompposition (Liu et e al., 2007; Yang ett al., 2001), esttimating depth by power specctrum

Mo odel parameter

Variaable name

Value

Sed diment conductiivity Sed diment heat production Bassement conducttivity Bassement heat production Mantle conductivity Mantle heat produuction Bassal temperature at the middle of the main dom Bassal temperature plan ne azimuth Bassal temperature plan ne inclination

Sediiment K Sedim ment HP

2.04 (W/(m·K)) 0..988 (µW/m3)

Baseement K Basem ment HP

2.31 (W/(m·K)) 1.35 1 (µW/m3)

Maantle K Manntle HP T0

3.28 (W/(m·K)) 1.28 1 (µW/m3) 1 075 (°C)

Azzimuth

28 (°)

Incllination

6.96e-4 (°)

The variable namee is the accessiible name in the t Guangdongg mod del input file to underworld, meeaning “Sedimeent K=3” at thee com mmand line set the value. Thhe basal temperrature plane iss effecctively 1 080 at a the SW cornner, 1 075 at the ES corner,, 1 07 74 at the NW coorner and 1 070 at NE corner.


85

Table 2 The simulated (red) and observational (black) geotherms of 14 of the boreholes of the model. In this simulation, the gradients and temperatures are largely co-incident, however, in certain boreholes there is a distinct offset. Some boreholes are affected by heat advection

0.5

3.5

4.0

0.0 0

Temperature Observed

-100

Depth (mAHD)

Guangzhou Conductivity (W/(m·K)) 1.0 1.5 2.0 2.5 3.0

Effective RL Simulated surface Conductivity

-200 -300

Depth (mAHD)

0.0 0

0.5

Lufeng Conductivity (W/(m·K)) 1.0 1.5 2.0 2.5 3.0

3.5

-10


-20


-30

4.0

-40 -50 -60

-400

-70 -80

280 290 300 310 320 330 340 350 Temperature (K) Huaxian Conductivity (W/(m·K)) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0


-100 -150 -200

Depth (mAHD)


-50

Depth (mAHD)

280 290 300 310 320 330 340 350 Temperature (K) Nanxiong Conductivity (W/(m·K)) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0

-100


-200


-300 -400 -500

-250 280 290 300 310 320 330 340 350 Temperature (K)

280 290 300 310 320 330 340 350 Temperature (K)

0.5

3.5

4.0

-100


-200


-300

0.0 0

Depth (mAHD)

Depth (mAHD)

0.0 0

Lechang Conductivity (W/(m·K)) 1.0 1.5 2.0 2.5 3.0


-200


-300

0.5

Puning Conductivity (W/(m·K)) 1.0 1.5 2.0 2.5 3.0

3.5

280 290 300 310 320 330 340 350 Temperature (K)

4.0

0.0 0


-100

-150

Depth (mAHD)

Depth (mAHD)


-200

4.0

-100

280 290 300 310 320 330 340 350 Temperature (K)

-50

3.5

-400

-400

0.0 0

0.5

Pingyuan Conductivity (W/(m·K)) 1.0 1.5 2.0 2.5 3.0

0.5

Renhua2 Conductivity (W/(m·K)) 1.0 1.5 2.0 2.5 3.0

3.5

4.0

-200


-400


-600 -800

-1 000 -1 200 280 290 300 310 320 330 340 350 Temperature (K)

280 290 300 310 320 330 340 350 Temperature (K)


86

Table 2 Continued

0.5

Depth (mAHD)

-50 -100

Qujiang Conductivity (W/(m·K)) 1.0 1.5 2.0 2.5 3.0

3.5

4.0

0.0 0


-40


Renhua1 Conductivity (W/(m·K)) 0.5 1.0 1.5 2.0 2.5 3.0

-100

0.0 0

4.0


-100 -150 -200

Depth (mAHD)

Depth (mAHD)

3.5

-250

280 290 300 310 320 330 340 350 Temperature (K) Shenzhen1 Conductivity (W/(m·K)) 0.5 1.0 1.5 2.0 2.5 3.0

3.5

-20


-40


-60

4.0

-80 -100 -120 -140

280 290 300 310 320 330 340 350 Temperature (K)

0.5

Shenzhen2 Conductivity (W/(m·K)) 1.0 1.5 2.0 2.5 3.0

3.5

4.0

0.0 0

Temperature Observed Effective RL Simulated surface Conductivity

-100

-150

-50

Depth (mAHD)

Depth (mAHD)

-80

-140


-50

-50

-60

4.0

-120 280 290 300 310 320 330 340 350 Temperature (K)

0.0 0

3.5


-200

-300

Shantou Conductivity (W/(m·K)) 1.0 1.5 2.0 2.5 3.0

-20

-150

0.0 0

0.5


Depth (mAHD)

0.0 0

-100

280 290 300 310 320 330 340 350 Temperature (K)

0.5

Wenyuan Conductivity (W/(m·K)) 1.0 1.5 2.0 2.5 3.0

3.5

4.0

Temperature Observed Effective RL Simulated surface Conductivity

-150

-200 -200

280 290 300 310 320 330 340 350 Temperature (K)

7 HEAT FLOW In its most basic configuration Underworld-GT implements the heat flow Equation (12) over a finite element computational domain supported by a particle based material domain. The heat conduction is driven through boundary conditions at the top and bottom of the domain taking the form of either a temperature or heat flux normal to the surface at each point of the surface’s mesh. Heat production and conduction is provided per particle in the domain and interpolated into the finite element domain to satisfy Equation (12) across the domain. The use of PETSc permits various direct and iterative solvers to be employed. These models used GMRES set to a relative tolerance of 1e-7. In the example case (13), the Underworld Picard iteration non-linear solver set to a tolerance of 1e-4, and

280 290 300 310 320 330 340 350 Temperature (K)

GMRES inner solves typically converges in a few iterations. The domain resolution is 185×136×87 elements in the eastings, northings and depth directions, calculated over 16 processors with 3 GB of memory per core. Each model, including the time to load the geometry, takes approximately 3 mins. The following is a particular realisation from an ensemble that satisfies some specific criteria generated by an inversion process to find likely material and basal conditions to meet observables. The specific criterion is that the simulated Curie depth reduced to a plane meets the observed Curie (described above) also reduced to a plane, with 500 km at each of the four corners of the domain. The plane comparison allows dimensional comparison between observation and simulated Curie. Given the simulated Curie is a surface of low curvature derived

Underworld-G GT Applied to Guangdong, G a Tool T to Explore the Geothermaal Potential of thhe Crust from the simpplicity of the geeological modell and boundary conditions this appproximation seeems to fit our need. In the exxperiments to saatisfy this criiterion, a feaature incremennt to Underworld-G GT occurred whhich permits thhe basal temperrature to be describeed as a bilinearr variation on thhe base. The suubsequent researchh question then became—whenn the Curie critterion is satisfied, arre bore hole geeotherms consiistent with obsservations? And dooes the basal tem mperature captture the fact thaat the Moho has a sttep in the domaain? These paraameters are not ideal but exemplify the progressionn of questioning. The results of o the total model with w these paraameters are shoown in Fig. 100 and Table 2. 8 CONCLUD DING REMAR RKS Shown heere is just one possible p param meter set from an ensemble of models that meet a collection off criteria includding a planar comparrison of Curie depths and tem mperature and gradig

877

c off observed geothherms. There are a some obser-ent comparisons vatio ons to be madee. For example this model hass not much dif-feren nce in conducttivity and heatt production between mantlee and basement. Is the 150º isothherm too low? Are the well-correelated boreholes those at distance from heat h advectionn proccesses? Is the surface heat fl flow reasonablee? The criteriaa must evolve to meeet the expectatioons from these questions. q To meet the crriteria, whetherr they be formaal or subjective,, g and changingg the modelling jourrney has requirred developing closeely related feaatures that are algorithms, daata input, dataa outp put and so on. Hence H the moddelling of Guan ngdong has fol-loweed the outlinedd approach. Furrther the tool Underworld-GT U T impllements this appproach, with ssome useful inp put and outputt alreaady implementted features foor regional scaale geothermall prob blems. The appproach is partticularly usefull when modell prov venance and sofftware reuse is iimportant.

Figure 10. Faacing west, usin ng the examplee parameters, shows the sim mulated Curie (gray) ( approximately over la apping the Cu-rie depth derrived through geophysical g m methods. Also shown is a 150 ºC isotherm. A visual comp parison of the simulated and d observation Curies C depth caan be made, raatifying inversiion results. ACKNOWLE EDGMENTS Underwoorld and Underrworld-GT has been supporteed by the Australiann Government infrastructure investment i in AuSA cope Simulation and Modellling and the NeCTAR Geologgy to Geophysics eR Research tool. Compute timee was supporteed by the National Computational C Infrastructure grant g “p67—Ennergy driven undersstanding of deeep geological resources” andd the NeCTAR Reesearch Cloudd grant “Undderworld_NeCT TAR_ Cloud_Flow” executed at thhe R@CMon, the t Research Cloud C node at Monassh University. Portions P of thiss paper appear in i the extended abstrracts to the Inteernational Workkshop of Deep Geothermal System ms, Wuhan 20112 and the Inteernational Workkshop on the Exploraation and Deveelopment of Geeothermal Resouurces in Guangzhouu, Guangdong 20013. REFERENCE ES CITED Balay, S., Eijkhout, V., Groopp, W. D., et al., 1997. Effiicient Managem ment of Parallellism in Object Oriented Numeerical Software Libraries. In: Arge, A E., Bruasset, A. M., Lanngtangen, H. P., eds., Moddern Software Tools in Scienntific Computinng. Birkhauser, Boston. 163–2202

ory, D., Sarvella, J. N., Rausschmayer, A., 2003. Scalingg Bato Step-Wise Refinement. In: P Proceedings of the 25th inter-national Confference on Software Engineering. Interna-tional Confereence on Softwaare Engineerin ng. IEEE Com-puter Society, Portland, May 3–10, 2003. 18 87–197 Capiitanio, F. A., Faaccenna, C., Zlootnik, S., et al., 2011. Subduc-tion Dynamicss and the Originn of Andean Oro ogeny and Boli-vian Orocline. Nature, 480: 833–86. doi:10.103 38/nature10596 Chop pra, P., Holgate, F., 2005. A GIS Analysis of o Temperaturee in the Australiian Crust. In: Prroceedings Worrld Geothermall Congress, Anttalya, April 24––29, 2005 Danis, C., 2014. Use U of Grounddwater Temperature Data inn Geothermal Exploration: E Thhe Example of Sydney Basin,, Journal, Australia. H Hydrogeology 22: 87–106.. doi:10.1007/s110040-013-1070-4 Danis, C., O’Neilll, C., Quenette, S., 2012. Is It Hot enoughh down there? Assessing G Geothermal Potential in thee Sydney-Gunneedah-Bowen B Basin. Australiaan Geothermall Energy Conferrence, Sydney. 14–16 Farr, T. G., Rosenn, P. A., Caro, E., et al., 200 07. The Shuttlee Radar Topogrraphy Mission. Rev. Geophys.., 45: RG2004.

88


doi:10.1029/2005RG000183 Gropp, W. D., Lusk, E., 1996. User’s Guide for mpich, a Portable Implementation of MPI. ANL/MCS-TM-ANL-96/6 Rev D, Mathematics and Computer Science Division Honza, E., Fujioka, K., 2004. Formation of Arcs and Back‐Arc Basins Inferred from the Tectonic Evolution of Southeast Asia since the Late Cretaceous. Tectonophysics, 384: 23– 53. doi:10.1016/j.tecto.2004.02.006 Kiczales, G., Lamping, J., Mendhekar, A., et al., 1997. AspectOriented Programming. In: Aks i̧ t, M., Matsuoka, S., eds., Proceedings of the European Conference on ObjectOriented Programming. Springer-Verlag, Berlin, Heidelberg, and New York. 1241: 220–242 Lafaurie, B., Nardone, C., Scardovelli, R., et al., 1994. Modelling, Merging and Fragmentation in Multiphase Flows with SURFER. Journal of Computing Physics, 113: 134–147 Li, C., van der Hilst, R. D., 2010. Structure of the Upper Mantle and Transition Zone beneath Southeast Asia from Traveltime Tomography. J. Geophys. Res., 115(B7). doi:10.1029/2009JB006882 Liu, T. Y., Wu, Z. C., Zhan, Y. L., et al., 2007. Wavelet MultiScale Decompostition of Magnetic Anomaly and Its Application in Searching for Deep-Buried Minerals in Crisis Mines: A Case Study from Daye Iron Mines. Earth Science—Journal of China University of Geosciences, 32: 31–41 (in Chinese with English Abstract) Ma, X. H., 1987. Outline of Chinese Lithospheric Dynamics. Geological Publishing House, Beijing (in Chinese) Mallet, J. L., 1992. GOcad: A Computer Aided Design Program for Geological Applications. In: Turner, K., ed., Three-Dimensional Modelling with Geoscientific Information Systems: Nato ASI Series C. Kluwer Acad. Publ., Dordrecht. 354: 123–141 Meixner, A. J., Holgate, F. L., 2009. The Cooper Basin Region 3D MapVersion 1: A Search For Hot Buried Granites. Geoscience Australia, Record 2009/15 Monaghan, J. J., 1992. Smoothed Particle Hydrodynamics. Annual Review of Astronomy and Astrophysics, 30: 534–574 Moresi, L., Dufour, F., Mulhaus, H. B., 2003. A Langrangian Integration Point Finite Element Method for Large Deformation Modeling of Viscoelastic Geomaterials. Journal of Computing Physics, 184: 476–497 Moresi, L., Quenette, S., Lemiale, V., et al., 2007. Computational Approaches to Studying Non-Linear Dynamics of the Crust and Mantle. Physics of the Earth and Planetary Interiors, 163: 69–82 Musson, A., Harrison, B., Gordon, K., et al., 2009. Thermal Thinking: Optimal Targeting for Australian Geothermal Explorers. In: Budd, A. R., Gurgenchi, H., eds., Proceedings of the 2009 Australian Geothermal Energy Conference. Geoscience Australia, Record 2009/35 Pavlis, N. K., Holmes, S. A., Kenyon, S. C., et al., 2008. An Earth Gravitational Model to Degree 2160: EGM 2008. General Assembly of the European Geosciences Union, Vienna Pohl, K., Böckle, G., van der Linden, F., 2005. Software Product Line Engineering—Foundations, Principles, and Techniques. Springer

Quenette, S., 2014. Software Evolution and Computational Performance—The Semantics of Software as an Inverse Problem: [Dissertation]. Monash University, Melbourne Quenette, S., Hodkinson, L., 2009. StgDomain—Scalable Parallel Domain Software Components for Particle-in-Cell Finite Element Methods. Concurrency and Computation: Practice and Experience, 22(12): 1593–1603 Quenette, S., Moresi, L., 2010. Models Based Experimentation: Numerical Modelling of 3D Basin Scale Architecture Heat & Fluid Flow. AGU Fall Meeting Abstracts, December 10, 2010, San Francisco Quenette, S., Moresi, L., Sunter, P. D., et al., 2007. Explaining StGermain: An Aspect Orientated Rnvironment for Building Extensible Computational Mechanics Modelling Software. HIPS Workshop, Parallel and Distributed Processing Symposium, Proceedings of the 19th IEEE International. 210. doi:10.1109/IPDPS.2007.370400 Somerville, M., Wyborn, D., Chopra, P. N., et al., 1994. Hot Dry Rock Feasibility Study. Energy Research and Development Corporation (ERDC) Report, 94/243, ERDC, Canberra Spector, A., Grant, F. S., 1970. Statistical Models for Interpreting Aeromagnetic Data. Geophysics, 35: 293–302 Stegman, D. R., Moresi, L., Turnbull, R., et al., 2008. Next Generation Visualization Framework for HighPerformance Computational Geodynamics. Visual Geosciences, 13(1): 71–84 Sulsky, D., Zhou, S. J., Schreyer, H. L., 1995. Application of a Particle-in-Cell Method to Solid Mechanics. Comput. Phys. Commun., 87: 236–252 USGS, 2004. Shuttle Radar Topography Mission. Global Land Cover Facility, University of Maryland, College Park, Maryland Wan, T., 2011. The Tectonics of China: Data, Maps and Evolution. Higher Education Press, Beijing. 173 Wand, M., Kiczales, G., Dutchyn, C., 2004. “A Semantics for Advice and Dynamic Join Points in Aspect-Oriented Programming,” ACM Trans. Programming Languages and Systems, 26(5): 890–910 Yang, W. C., Shi, Z. Q., Hou, Z. Z., et al., 2001. Discrete Wavelet Transform for Multiple Decompostion of Gravity Anomalies. Chinese Journal of Geophysics, 44: 534–541 Zhang, B. T., Wu, J. Q., Ling, H. F., et al., 2007. Estimate of Influence of U-Th-K Radiogenic Heat on Cooling Process of Granitic Melt and Its Geological Implications. Science in China Series D: Earth Sciences, 37: 155–159 Zhang, Z., Xu, T., Zhao, B., et al., 2013. Systematic Variations in Seismic Velocity and Reflection in the Crust of Cathaysia: New Constraints on Intraplate Orogeny in the South China Continent. Gondwana Research, 24(3): 902–917 Zhou, L., Xie, J., Shen, W., et al., 2012. The Structure of the Crust and Uppermost Mantle beneath South China from Ambient Noise and Earthquake Tomography. Geophysical Journal International, 189(3): 1565–1583 Zhou, X., Sun, T., Shen, W., et al., 2006. Petrogenesis of Mesozoic Granitoids and Volcanic Rocks in South China: A Response to Tectonic Evolution. Episodes, 29(1): 26