Mathematical Geosciences: Local Singularity Analysis of Nonlinear Earth Processes and Extreme Geo-Events

In the first part of the chapter, the status of the discipline of mathematical geosciences (MG) is reviewed and a new definition of MG as an interdisciplinary field of science is suggested. Similar to other disciplines such as geochemistry and geophysics, mathematical geosciences or geomathematics is the science of studying mathematical properties and processes of the Earth (and other planets) with prediction of its resources and changing environments. In the second part of the chapter, original research results are presented. The new concepts of fractal density and local singularity are introduced. In the context of fractal density and singularity a new power-law model is proposed to associate differential stress with depth increments at the phase transition zone in the Earth’s lithosphere. A case study is utilized to demonstrate the application of local singularity analysis for modeling the clustering frequency—depth distribution of earthquakes from the Pacific subduction zones. Datasets of earthquakes with magnitudes of at least 3 were selected from the Ring of Fire, subduction zones of Pacific plates. The results show that datasets from the Pacific subduction zones except from northeastern zones depict a profound frequency—depth cluster around the Moho. Further it is demonstrated that the clusters of frequency—depth distributions of earthquakes in the colder and older southwestern boundaries of the Pacific plates generally depict stronger singularity than those obtained from the earthquakes in their hotter and younger eastern boundaries.

Congress (IGC) in Prague. It had grown from mathematical geology to mathematical geosciences by the time its name was changed at the 32th IGC held in Oslo in 2008. Not only has the subject been accepted widely within the geoscience community but the association has also been recognized for its reputation and significant influence on the earth sciences in general. IAMG has affiliations with several major geoscience organizations including the International Union of Geological Sciences (IUGS), International Statistical Institute (ISI), and the International Union of Geodesy and Geophysics (IUGG). Diverse earth science topics have been published in IAMG conference proceedings and the IAMG journals (Mathematical Geosciences, Computers & Geosciences and Natural Resources Research). However, we have to realize that as a relatively young discipline, MG still has not been very widely accepted and is often ignored by main stream geoscientists. While several definitions and terminologies were proposed to describe mathematical geology, there have been few attempts to define mathematical geosciences. For example, mathematical geosciences have often simply been referred as applications of mathematical and statistical methods for the analysis of geological (earth science) data and the development of quantitative predictive models (Howarth 2017). The mission of the IAMG as shown on the IAMG website was defined as promoting the development and application of mathematics, statistics and informatics in the geosciences. Whether MG should be defined as a formal discipline of science or simply as applications of mathematics in the geosciences is a fundamental question with critical impact on the development of the subject. In this chapter, I will review the status of the discipline and suggests a new definition for MG followed by examples to demonstrate what contributions of MG have been made to earth science and what the current developments in the field are. For the first part I will elaborate on MG on the basis of literature review and for the second on my own research in nonlinear MG as an example of a new field of MG.

What Is Mathematical Geosciences or Geomathematics?
One of the original definitions of mathematical geology was given by Vistelius (1962) and used in the name of the association: International Association for Mathematical Geology (IAMG) when it was first established in 1968. Geostatistics is one of the successful fields of IAMG, which originally was developed by MG scientists within the IAMG community. It has been used not only in the geosciences but later in many other fields of science as well. Geostatistics focuses on application of statistical methods in the earth sciences (e.g. Merriam 1970;McCammon 1975a, b) and still appears to be used by many in that sense. The term geomathematics was also used by several authors including Agterberg (1974) who used the term as the title of his two books (Agterberg 1974(Agterberg , 2014. After the name of the association was changed from mathematical geology to mathematical geosciences in 2008, the term mathematical geosciences more often appears in the literature of the IAMG and also in the titles of conferences, as well as in the name of its journal Mathematical Geosciences. When the author of the current chapter served as president of IAMG (2012)(2013)(2014)(2015)(2016), dedication to IAMG was given by promoting the discipline of mathematical geosciences. Several notes on this were published in the President's Forums in IAMG newsletters . The distinction between mathematical geology and mathematical geosciences is not simply in terminology but also in the scope of the discipline. While mathematical geology refers to a branch of geology, mathematical geosciences must be a subdiscipline of the geosciences which includes geology as one of its subfields. Other relevant subjects covered in the geosciences include but are not limited to geochemistry, geophysics, geobiology, and hydrology. Mathematical geosciences should be a discipline parallel to other subdisciplines in the geosciences such as geochemistry, geophysics and geobiology rather than a branch of geology. In the author's personal view this distinction is critical for the development of the discipline. Under the concept of mathematical geology, the subject is limited to the application of mathematics in geology but as mathematical geosciences just like geochemistry and geophysics, it serves the entire earth science. So, what should be the definition of mathematical geosciences or geomathematics and what are the roles mathematical geosciences should play in the family of geosciences? Here I will briefly elaborate on these questions and introduce several major contributions of MG to earth science. In order to provide a proper definition of mathematical geosciences, we should look at the definitions of other relevant disciplines such as geochemistry, geophysics and geobiology: • Geophysics as a science of "the study of the earth's physical properties and of the physical processes acting upon, above, and within the earth." (Collins English Dictionary) • Geochemistry as a science that deals with the chemical composition of and chemical changes in the solid matter of the earth or a celestial body (Unabridged dictionary). • Biogeosciences as an interdisciplinary field of study integrating geoscience and biological science: the study of the interaction of biological and geological processes (Unabridged dictionary).
The definitions of the preceding relevant disciplines share the common concept of an interdisciplinary geoscience field. A similar definition was proposed by the author in 2014 with consultation of the IAMG Executive Committee Members and published in the President's Forum of IAMG newsletter (Issue No. 79).
• As an interdisciplinary field merging mathematics, computer science and geosciences, MG is the science of studying mathematical properties and processes of the Earth (and other planets) with prediction and assessment of its resources and environments The ultimate question arising from this definition is what are the mathematical properties and processes of the Earth, with prediction and assessment of its resources and environments which have to be dealt with by mathematical geoscience for integration with other geoscience subdisciplines. Similar to other interdisciplinary fields including geochemistry and geophysics, mathematical subjects such as geometry, calculus, functional analysis, morphology, probability and mathematical statistics provide essential theory and methods for quantitative study of the Earth ranging from geometry and dynamics of the Earth, uncertainties of measurements, and observations for the prediction of Earth events.

What Contributions Has MG Made to the Geosciences?
There are many examples demonstrating that MG has made indispensable contributions to the geosciences. For example, the mathematical model of the Earth's shape (e.g. Clark ellipsoid, and Hayford ellipsoid) which serves as the foundation of geodesy, navigation systems (e.g. GPS), remote sensing technology (RS) and geographical information systems (GIS), and the fast growing field of geomatics; the mathematical model of mantle convection and models for plate motions (McKenzie and Parker 1967) serve as foundation of plate tectonics, the most notable development of earth science in the last century; mathematical symmetry and symmetry operations as principles of crystallography and optical mineralogy (e.g. in 1830, Hessel proved the existence of the 32 groups of crystal symmetry) which constitute a foundation of solid earth science; the mathematical topological model as foundation of geographical information systems (e.g. as basis of spatial data modeling in ArcGIS), one of the most useful technologies in geoscience; mathematical and statistical theories providing foundations for describing the spatial distribution and correlation of elements, uncertainty and error bars in geochemistry including isotope geochemistry and geochronology as are also used for the geological time scale; mathematical modeling and uncertainty of prediction of climate change, a pressing issue of the geosciences; probability theory and stochastic models for prediction of energy and mineral resources, highly demanded by many nations for economic and societal development; geo-complexity theory such as fractals, multifractals, chaos and self-organized criticality for modeling and predicting singular events and extreme phenomenal issues; and information extraction (big data mining, machine learning, geo-intelligence) in the geosciences, just to name a few. As the International Association for Mathematical Geosciences, IAMG has earned its reputation by promoting and fostering its members to make contributions to science. Original and significant studies have been published in IAMG journals, books and conference proceedings. However, a large amount of work is documented elsewhere in publications which cover almost every mathematical subject and aspects of geosciences ranging from statistical data analysis, geometrical modeling, dynamics and processes simulation, to prediction and assessment of Earth system. MG theories and methods have been applied not only in tackling conventional solid earth issues such as assessment of mineral and energy resources, but also in other fields including hydrology, climate change, water resources, alternative energy resources and environmental issues. While the importance of MG in the geosciences has been increasingly demonstrated, the discipline of MG has not yet been fully recognized and, to some extent, buried in oblivion. There is hardly any hiring of highly qualified personal (HQP) in academic institutions or industry with as job title Mathematical Geoscientist or Geomathematician. As a matter of fact, most of our IAMG members are employed with job titles such as geologists, geophysicists, geochemists, geodesists, computer scientists, mathematicians and geoinformatical specialists instead of MG or GM. University students who are talented in mathematics and geosciences wanting to pursue mathematical geoscience have to enroll in geophysics or other fields simply because MG does not exist as such in university programs, at least in most of the programs in developed nations. There are very few interdisciplinary university programs except actuarial science, mathematical physics and mathematics for business, which have mathematics as integral part of their subject. A common misconception is that learning mathematics either can only result in kinds of two jobs: pure mathematician or mathematics teacher, or as a prerequisite for other careers in engineering, science or business. This might be one of the reasons there are not so many students wanting to pursue mathematics related subjects in their choice of career. Thus, MG faces significant challenges when promoting MG as a discipline and for facilitating training and education of future generations. This presents the bottleneck for the IAMG to grow further and to become a more successful and influential association.
The International Year of Mathematics of Planet Earth (MPE) celebrated in 2013 generated a much needed publicity of mathematics in geoscience. Mathematical courses are offered in all schools from primary to high school to university. Earth science is also a common choice of topic in essays by students. Integration of math and earth subjects must provide proper and interesting topics for students' math or science projects. The mathematical and geoinformatical techniques learned by students early on are already powerful tools for exploring the Earth. An excellent example is the work headlined in the media with publication by a high school student Alice R. Zhai who has analyzed 73 tropical cyclones that made landfall in US and used multivariate regression to examine the dependence of hurricane economic loss on maximum wind speed and storm size. This study (Zhai and Jiang 2014) not only proposes a new model by which hurricane damage might be predicted but also provides new evidence showing the area-density power law property of extreme events which, as is to be introduced in the remainder of this chapter, has deep origins in nonlinear dynamics.
The development of modern information technology enables everyone to easily retrieve big data to support their studies via internet and web services in a cloud environment. To access and process huge amounts of data is no longer only for paid professionals. More and more specialized software packages and multi-media teaching and training materials or online courses available in the public domain with Twitter, Facebook and You Tube, provide new ways for self-learning. Online communication, discussion and consultation through the internet in and out of the classroom have become common for students. It should encourage middle school, high school and university students to develop their curiosity in, passion for, and dedication to mathematical geosciences.

Frontiers of Earth Science and Opportunities of MG
IAMG has been rapidly expanding its scope from traditional geostatistics or statistical geology to more comprehensive interdisciplinary sciences for mathematically studying properties and processes of the Earth with prediction and assessment of its resources and environments. What are the current trends of MG and how are they associated with the Earth Science frontiers? It is impossible to create an accurate list of frontiers for MG. Of course, there exist several previous publications by IAMG members that have discussed past, current and future trends for the IAMG (Agterberg 2003). Here I will just share some thoughts based on my personal observations of several recent events and activities. Several international organizations have developed and published white papers illustrating prospective review on trends of scientific research within their organizations and strategic plans for the next 5-10 years; for example, the International Council for Science Union (ICSU) published its strategic research agenda for Future Earth 2025 Vision (http://www.futureearth.org/ sites/default/files/future-earth_10-year-vision_web.pdf); the International Union of Geological Sciences (IUGS) is jointly with UNESCO offering the International Geological Correlation Program (IGCP) in addition to various other big science programs and new initiatives such as the Resourcing Future Generations (RFG), an international collaborative program (http://iugs.org/uploads/RFG.pdf); the US National Science Foundation (NSF) has published a strategic plan for 2014-2018 (https://www.nsf.gov/publications/pub_summ.jsp?ods_key=nsf14043); the American Geophysics Union (AGU) produced a scientific trends report (https://about.agu. org/trends-earth-space-science/); the American Natural Science Foundation published its strategic plan for tectonics (https://docgo.net/national-science-foundationnsf-strategic-plan-fy-2006-2011-nsf-06-48); a white paper resulting from NSF sponsored workshops on "mathematics in geosciences" was published by a group of geoscientists in 2012 (https://cpb-us-e1.wpmucdn.com/sites.northwestern.edu/dist/ 8/1676/files/2017/10/agenda-xwphux.pdf), just to name a few. Relevant publications resulting from international conferences such as the International Geological Congress (IGCs), AGU, EGU, GSA as well as special articles in several journals such as Nature and Science have also been concerned with these issues. The following summary of key topics can be extracted from the preceding sources of information to reflect current trends and frontiers of the earth sciences. These key topics include but are not limited to data science, data analysis, big data and geo-intelligence, computation, inter-/multi-/cross-/transdisciplinary science, integrated models, uncertainty relative to observations and predictions, properties and dynamics of the planet, climate change, disruptive processes such as earthquakes and storms, and special studies of the Arctic, Antarctic and Tibet Plateau. The fundamental issues are for understanding Earth and environmental systems and their interactions with human activities, and for developing reliable monitoring systems, models, and information technologies for predictions and early warnings of large-scale and rapid change. The current challenges facing earth scientists are understanding and modeling the geo-complexity of the Earth and environmental systems with their interactions, chaotic nature and predictability of geo-processes, Earth singularity and human mitigation and adaptation to extreme events, plus observation and monitoring multiple-scale mixing nonlinear processes. Although most organizations neither recognize nor explicitly mention this, the majority of these frontiers are fundamentally related to MG. A long period of incremental advances of new mathematical theories and models in conjunction with modern technologies for solving these earth science problems may lead to creative leaps of innovation. MG has huge challenges and responsibilities facing the earth science frontiers. MG scientists are indeed at the frontier of earth science tackling fundamental problems of the Earth as can be evidenced by the recent advancements reflected in the topics of plenary presentations at IAMG conferences and in the best papers published in IAMG journals; for example, on multi-point geostatistics-a new field of spatial-temporal modeling (Mariethoz and Caers 2014); compositional data analysis-a new way to explore the composites of the Earth (Pawlowsky-Glahn et al. 2015); singularity analysis and singularity physics-new theory and methods of studying geodynamics and geo-complexity (Cheng 2007(Cheng , 2017aAgterberg 2017); big data visual analytics for exploratory data analysis; semantic web technology for geoinformation; uncertainty in ecosystem mapping by remote sensing; integrating structural geological data into inverse modeling frameworks; stationary and isotropic vector random fields on spheres; and mathematical morphology modeling, just to name a few.

Fractal Density and Singularity Analysis of Nonlinear Geo-Processes and Extreme Geo-Events
For the past several decades nonlinear theory and geocomplexity marked an era of new geoscience that deals with nonlinear processes and extreme phenomena which occurred in the evolution of earth systems. Irregular geometry was not popularized in the past until the term "fractal" was coined by Mandelbrot in the 1970s. Fractal geometry rapidly became a new field of mathematics dealing with roughness and irregularity of geometries. For example, fractals have been used for modeling complex and self-similar patterns generated by nonlinear processes (Mandelbrot 1972;Feder 1988). The concept of fractals and fractal dimension was further extended to multifractals involving self-similar measures defined on support which can be fractal itself (Mandelbrot 1972;Meakin 1987;Schertzer and Lovejoy 1987). Multifractal measures have been further extended to fractal density in local singularity analysis (Cheng 1999a(Cheng , 2001. In the following sections the concept of fractal density will be introduced and followed by discussion and application of new methods for fractal differential operation and fractal integration (Cheng 2017a).

Fractal Density
Since the principle of density was discovered by the Greek scientist Archimedes approximately 2000 years ago, the well-known physical concept of density has become a fundamental property of mass or energy with a variety of applications. The density, or volumetric mass density, of a substance is its mass per unit volume. Density thus is a scale-independent property of material or energy treated as representing a fundamental physical parameter and variable in many physical models with applications in nearly all fields of study, ranging from physics to engineering, economics and the social sciences. Density often is characterized as unit of mass over volume (e.g., g/cm 3 , kg/m 3 ) or energy over volume (J/cm 3 , w/L 3 ). For example, the density of pure gold is 19.32 g/cm 3 , which is approximately 19 times as much as for an equal volume of water. The density of quartz is 2.65 g/cm 3 , which is much less than the density of gold. Therefore, the density of gold-mineralized quartz veins in hydrothermal mineral deposits is variable depending upon the concentration and distribution of gold in the quartz veins. Similarly, continental crust, which consists mostly of granitic rock, has a density of about 2.7 g/cm 3 and the Earth's mantle of ultramafic rock has a density of about 3.3 g/cm 3 . The density of seawater varies with temperature and salinity of the water. Although the density of seawater varies at different points in the ocean, a good estimate of its density at the ocean's surface is 1025 kg/m 3 or 1.025 g/cm 3 . Density of air is a temperature and pressure dependent parameter. For given temperature and pressure the density of air is independent of the volume of air. For a pure substance the density is independent of the volume of substance. However, for a heterogeneous substance density usually assumes different values depending upon purity and packaging. For example, rocks consisting of minerals with different densities have variable densities depending upon the proportions of the minerals. For a quartz vein with pure SiO 2 the density of the vein should be equal to the density of quartz, 2.65 g/cm 3 . However, if the quartz vein involves gold mineralization, then the density of the quartz will be different from that of pure quartz relating to how the gold is distributed in the vein. At a location of higher concentration where a cluster of gold occurs in the quartz vein, the density of the vein is higher than that of pure quartz. From a fractal point of view, the structure of these types of gold distribution can be very irregular and then has to be described by using a non-integer or fractal dimension. Accordingly, the value of "volume" of the substance is lost. Instead the size of fractal is measurable only if it is measured in fractal dimensional space or as Hausdorff measure (Cheng 2017a). This means the ratio of mass over volume does not converge; and the density does not exist according to the ordinary density definition. In the following section it will be demonstrated that the concept of ordinary density of substance is only valid for substances with regular or ordinary structure. For substances packaged in a fractal manner, a new form of density is needed and the concept of ordinary density has to be generalized to a new form capable for quantifying the density of complex objects. It will also be demonstrated that the end products for many types of singular processes possess fractal mass density or energy density. The concepts of fractal density and local singularity analysis have been utilized in several dynamic models involving extreme processes (Cheng 2012(Cheng , 2016(Cheng , 2017bCheng and Agterberg 2009;Cheng and Sun 2017).

Density-Scale Power-Law Model and Singularity
According to the concept of ordinary density, the mass density of an object (ρ) can be calculated by the following equation: where m(v) represents the mass contained in a volume (v) and ρ is the average density of an object. If the object is homogenous then the density calculated in Eq. (10.1) becomes independent of volume. The unit of the density is determined by the ratio of the mass and volume; for example, g/cm 3 . However, if the object has heterogeneous properties, the density may vary from place to place and the average density in Eq. (10.1) varies with different size of v, then a localized density must be calculated using the derivative of the mass over volume: The density in Eq. (10.2) exists only if the limit converges when the volume becomes infinitesimal. If the limit does not converge, then the density doesn't exist. As a generalization of Eq. (10.2), the following new Eq. (10.3) was introduced (Cheng 1999b(Cheng , 2001 in which there exists a parameter α (with positive value) so that the limit converges: .

ð10:3Þ
The value of ρ α can be considered as a generalized density because the ordinary density defined in Eq. (10.2) becomes a special case of Eq. (10.3) when α = 3, the normal dimension of volume. This new density was named fractal density since it is defined as mass or energy per unit of "fractal set" (Cheng 1999b(Cheng , 2001. The fractal density defined in Eq. (10.3) has as unit the ratio of mass to a fractal set of α dimensions; for example, g/cm α or kg/m α . Similarly, the units of fractal energy density can be J/cm α or w/L α . Combining Eqs. (10.2) and (10.3) yields the following relationship between ordinary density and fractal density: ð10:4Þ The notation of fractal density used in Eqs. (10.3) and (10.4) can be replaced by the following general model associating the fractal density with the ratio of mass and scale (ε-linear size of an E-dimensional set): ð10:5Þ This power-law relation between the ordinary density and scale is determined by two parameters: the fractal density ρ α which is independent of scale and the exponent-singularity index α (fractal dimension), or Δα = E − α; the latter is also known as the co-dimension of fractal density. The singularity index (Δα) measures the deviation of the fractal dimension from the dimension of normal density. These two parameters (ρ α and Δα) can be estimated from observed data by measuring the intercept and slope of a straight line on the log-log plot of m against ε (Cheng 1999b(Cheng , 2007.

Multifractal Density
If fractals refer to geometry with irregular shapes and self-similar geometrical properties, multifractals refer to self-similar measures defined on support which can be fractal (Mandelbrot 1983). Multifractals are defined as spatially intertwined fractals with variable fractal dimensions (e.g., Mandelbrot 1972;Cheng 1997). According to the distribution of measures (similar to the mother functions of sets) the support can be grouped into subsets which can be fractal with specific fractal dimension. Accordingly, there are two types of multifractal measures: continuous and discrete multifractals, the former refers to multifractals corresponding to an infinite number of intertwined fractals with continuous fractal dimension spectrum, whereas the latter refers to the limit number of intertwined fractals with discrete fractal dimensions (Cheng 1997). Multifractal measures are self-similar measures with multiple scale singularities which can be characterized by the Hőlder exponent (Mandelbrot 1989). In the multifractal paradigm the measure defined on a support can be expressed as where μ(ε) represents the measure defined on a set of linear scale ε, ∝ stands for ''proportional to'' when cell size ε approaches to zero, and α is the singularity index also known the Hőlder exponent (Mandelbrot 1989). This power law exists usually in a statistical sense and is represented as expectation <>. According to the distribution of α values, the entire support can be classified into subsets or fractals each with different singularity and accordingly different fractal dimensions. This is why it has been termed "multifractal". The distribution of singularity α in the mapped area can be described by the fractal dimension spectrum function f(α). The values of singularity and multifractal spectra can be estimated by several methods including box-counting and gliding-box based moment methods, and the wavelet method (Cheng 1999b). Singularity property has been commonly observed in geochemical and geophysical quantities (Cheng et al. 1994;Cheng 1999bCheng , 2007. Since the common moment-based multifractal models are implemented according to partition functions of measures with additive property, most literature about multifractals focuses on the power law relations of multifractal measures and self-similarity of multifractal measures and few have neither emphasized the physical meaning nor the property of density of the multifractal measure. A density-area fractal model was proposed (Cheng et al. 1994) to associate the concentration with area of multifractal measure as where the area (A) is a function of element concentration above the threshold C. The model has also been applied to characterize other types of "concentration" such as density of faults per area , density of mineral deposits per area , stream density per drainage area , and digital number of remote sensing images (Cheng and Li 2002), just to name a few. Further utilizing the idea of C-A model locally, the following power law relation was introduced to associate the density of multifractal measures with scale (Cheng 1999b) where E is the Euclidean dimension of the support (e.g., E = 1 for line, 2 for area and 3 for volume), x indicates the location, and c(x) and α(x) are constants with respect to scale ε but varying with location. The values of α(x) and c(x) can be estimated from the values ρ ε, x ð Þcalculated for different sizes ε around the location x by means of least squares using log-log paper. Both values can be mapped for visualization and interpretation. For convenience without loss of generality, in the rest of the paper the notation of x will be dropped from the formulation and the equation is assumed to hold locally. The singularity index α and constant c have the following properties (Cheng 1999b): if α = E, then ρ ε ð Þ = constant, independent of vicinity (scale) size ε; if α > E then ρ ε ð Þ is a decreasing function of ε which implies the convex property of ρ ε ð Þ; and α < E then ρ ε ð Þ is an increasing function of ε which implies the concave property of ρ ε ð Þ. Thus, the ordinary density obeys a power-law relationship with scale which has the following properties (Cheng 1999b(Cheng , 2007: ð10:9Þ In accordance with these properties, ordinary density becomes volume dependent when α ≠ E and it tends to either zero or infinity when the scale ε becomes infinitesimal. The constant c in Eq. (10.8) can be expressed in the following form: The constant c indeed is a convergent value of the ratio of measure (μ) over scale (ε) with fractal dimension. This quantity is usually termed scaling factor but it can be termed as a fractal density or Hausdorff density in analogy to the mass density which corresponds to ratio of measure over ordinary geometry with integer dimension (Cheng 2015). Therefore, while a unit of ordinary density is g/m E , the unit of fractal density becomes g/m α .

Fractal Density Structure and Clustering Distribution
The terminology of fractal density has been explained in several papers with different emphases, but the meanings of the concepts used are variable. For example, the term "fractal density" has been used to refer the number of fractals per area (Hou and Wu 1989) which does not mean the same as the concept introduced in the current paper. Tatekawa and Maeda (2001) analyzed time evolution of fractal density perturbations in the Einstein-de Sitter universe, in which the emphasis is on how the perturbation evolves and what kind of nonlinear structure will come out. Similarly, Federrath et al. (2009) has used fractal density structure in supersonic isothermal turbulence when referring to density structure. Gromov et al. (2001) used fractal density to describe fractal galaxy distribution. Carpinteri et al. (2009) used the term to describe the mean fractal density of microcrack barycenters. Pope and Mackenzie (1988) introduced the concept of fractal density for describing the morphology of fractal growth model in the evolution of gels from solution. They define the fractal density ρ which follows the relation where D is the fractal dimension of fractal growth, the F is the relative fractal density at radius r (r ≥ r 0 ), with r 0 and ρ 0 being the core radius and core density, respectively. The core acts mathematically as a reference point for calculating the decrease in density as the fractal increases in size. A similar clustering fractal growth density function was used to describe tumor growth in fractal space-time with temporal density (Paramanathan and Uthayakumar 2011).
From the preceding publications we can see that in earlier studies by other authors the term of fractal density was introduced mainly for description of morphology and patterns of fractals and fractal growth modeling. The current research introduces the fractal density as a generalization of ordinary density of substance or energy to represent a fundamental new parameter or variable involved in dynamic systems.

Fractal Integral and Fractal Differential Operations of Nonlinear Functions
As mentioned in Eq. (10.2) for heterogenetic matter or substances, the derivative of mass over scale can be used for defining localized density of substance. Accordingly, the mass or volume of a heterogenetic substance can be calculated using integration. Obviously, integration and differentiation are two fundamental operations in calculus and used for many mathematical and physical subjects. The traditional integral and differential operations are defined on the basis of additive property of Lebesgue measure. When the measure no longer possesses additive property, then the classical integral and differential may not exist. Therefore, the ordinary integral and differential operations are not applicable to fractal density with singularity. The author has proposed the following fractal integral and differential (Cheng 2017a) where Δf(x) and Δx represent the increments of a function f(x) for an increment of x. The convergence of the limit in Eq. (10.12) can be defined as the α-fractal derivative of the function f(x). Similarly, we can define the fractal integral of the function f(x) as follows where f(x i ) is the magnitude of the function f(x) over the small range [x i , x i + Δx].
If the limit of Eq. (10.13) converges, then it can be named the α-fractal integral of the function f(x). It must be kept in mind that the fractal derivative defined in this paper is different from the fractional derivative (fractional order) known in the literature as f (v) (x), where v can be a non-integer order. The fractional derivative assumes that the normal integer order derivative f (n) (x) does exist. The fractal derivative is based on fractal dimension of the measure whereas the fractional derivative is based on fractional order of derivative defined on normal measure. As an example, let us take a power-law function to demonstrate the fractal derivative. Assume a power law function, f(x) = c(x − x 0 ) b , with ordinary derivative of the function f ′ (x) = cb(x − x 0 ) b−1 , which does not exist at x = x 0 if 0 < b < 1. The integral of the function then is R f ðx) dx = c ̸ ðb + 1Þðx − x 0 Þ b + 1 , which does not converge if b < −1 at x = x 0 . According to Eq. (10.12), the fractal derivative at x = x 0 exists and f α ′( A new concept of Hausdorff derivative underlying the Hausdorff dimension of metric space/time was proposed by Chen (2006) who introduced the systematic mathematical operation of Hausdorff derivative with applications to derive a linear anomalous transport-diffusion equation underlying an anomalous diffusion process. The Hausdorff derivative operation proposed by Chen (2006) is expressed as follows This formalism was termed the Hausdorff derivative of a function f(x) with respect to fractal measure x α .
It has to be pointed out that the fractal derivation defined in Eq. (10.12) is different from that defined in Eq. (10.14) considering that, in general, if x 0 ≠ 0, then The two sides in Eq. (10.15) become equal only if x 0 = 0. Otherwise, according to Taylor expansion, we can obtain Δx α = x α − x α 0 = αx α − 1 0 Δx + o Δx ð Þ, so substitution into Eq. (10.14) gives

Earth Dynamic Processes and Extreme Events
In the remainder of this chapter I demonstrate that fractal density (Δα ≠ 0) characterizes anomalous mass accumulation or energy release caused by extreme geo-processes, which occurred in the Earth's lithosphere originated from cascade earth dynamics (plumes, mantle convection and plate tectonics) and self-organized criticality involved in phase transitions (avalanches of slab breakoffs, faults, and volcanic eruptions). Mantle convection at high Rayleigh number generates thermal plumes episodically which upon arrival in the crust could cause major flood basalt events, igneous provinces as well as spreading of continents and mid ocean ridges (Richards et al. 1989;White and McKenzie 1989). On a larger scale, Wilson cycles (Wilson 1966) corresponding to the periodic fragmentation and reformation of supercontinents could be linked to temporal variability in plate tectonics. Numerous studies have revealed that mantle convections can induce exchange of mass between upper and lower mantle across the endothermic phase transition zone at about 660 km. The cold downwelling material penetrates into the lower layer and, simultaneously, the hot upwelling fluid is pushed into the upper layer. The exchange of mass between the upper and lower mantle layers can occur in short bursts (often described with superlatives such as "catastrophic", overturn, "avalanche" subduction, or "superplumes") (Zhong and Gurnis 1994). The quick injection of lower mantle hot fluid into the upper mantle can cause not only mantle heterogeneity but also anomalous thermal distribution near the surface (Le Bars and Davaille 2004). This has been considered to be the first order cause of vigorous magmatism. Deep subductions of continental crust into the deep earth interior and rebounded back to the surface of the Earth have been ascertained by the discoveries of regional metamorphic coesite (Chopin 1984;Smith 1984), and subsequently by unusual ultrahigh pressure (UHP) terranes (Hacker and Gerya 2013).
Within the lithosphere there are various types of "catastrophic" events occurring during plate subduction. Formation of magmatic arc can be caused by subduction in which the subducting or subducted oceanic crust material releases volatiles (e.g. H 2 O and CO 2 ) which cause partial melting of the mantle and form magma at depth under the overriding plate. Earthquakes occur at certain depths at the edges of three types of plate boundaries: convergent (subductions and collisions), divergent, and transformative.

Phase Transition
From mathematical and physical points of view, the mechanisms that have been proved to exist correspond to the generation of power-law distributions including but not limited to phase transition (PT), self-organized criticality (SOC) and multiplicative cascade processes (MCP) (Newman 2005;Lovejoy et al. 2009). I will elaborate on each of these mechanisms in relation to mantle convections, plumes and lithosphere rheology induced tectonic events. The phase of a thermodynamic system and the state of matter in a normal system have uniform physical properties. Common phases include liquid phase, solid phase and vapor phase of chemical components which exist under certain pressure and temperature (P-T) conditions. Materials in different phases have their distinct properties such as liquid usually having higher density and smaller specific volume in comparison with gas. However, in phase transition conditions, multiple phases coexist within the same system such as liquid and vapor in magma and hydrothermal systems under proper P-T conditions. At a critical condition (critical point on phase diagram) liquid and vapor become indistinguishable and beyond this point the fluid and gas become so-called supercritical fluid, representing a special phase of matter which can effuse through solids like a gas, and dissolve materials like a liquid (McMillan and Stanley 2010). The critical point for water occurs at temperature (374°C) and pressure (22 MPa). It has been found that the critical point is so peculiar that close to it, small changes in pressure or temperature result in large changes in density and other density related properties such as viscosity, relative permittivity, heat capacity and solubility. The special critical point phenomena can be expressed by the following empirical power law functions Levelt Sengers 1968, 1986): where Δρ, ΔP, and ΔT represent the changes of density, pressure and temperature, respectively along the coexistence curve. These power-law relations hold for small changes of temperature or pressure from the condition at the critical point of the system. Although the two functions of Eq. (10.17) show continuity at zero increment with Δρ = 0, ΔP = 0, and ΔT = 0, the first order derivatives of density versus either temperature or pressure (change rate of density difference) do not exist or show singularity at ΔP = 0 and ΔT = 0 as shown in the following forms These properties describe the phenomena of property change such as fractal density (density jump) at the phase transition zone. In addition, the ratio of increments of temperature and pressure depict power-law relations ΔP ΔT = cΔP − 1 ̸ 3 . Such power-law relation implies that the Clapeyron slope could become infinity or a singularity when approaching the coexistence curve. Clapeyron slope and density jump are critical parameters in numerical simulation of mantle convection; for example, Korenaga (2004) developed a numerical model to simulate mantle mixing and continental breakup magmatism by assigning a Clapeyron slope of −2 MPa/K and a density jump of 10% for the endothermic phase transition at 660 km depth. The episodicity of convection induced by the endothermic phase changes strongly depends on plate length, rheology, and Clapeyron slope (Zhong and Gurnis 1994). Ogawa and Yanagisawa (2014) have developed models with small Clapeyron slope −0.2 to −1 MPa/K for simulating convections from punctuated layered convection to whole-mantle convection in modeling mantle evolution on Venus due to magmatism and phase transitions. Their models indicate that the earlier stage layered mantle convection is punctuated by repeated bursts of hot material from the deep mantle to the surface. Other phenomena of phase transition may occur at the boundary of deeply subducted slabs. Due to subduction of oceanic lithosphere underneath the continental lithosphere, solid phase lithosphere can be partially melted to facilitate formation of magma. During the progress of subduction, H 2 O and other volatile components contained in the rocks are progressively released from the slab at different depths. Fluids or melts released at greater depths will be in supercritical fluid phase which hydrates the mantle and causes partial mantle melting. This eventually leads to deeply rooted magma which provides the source for magmatic and volcanic arcs located above the subduction zones. Partial melting in lower crust and mantle also causes strain rate change of the lithosphere which facilitates formation of intermediate and deep earthquakes (Dimanov et al. 2000). The processes of fluid release and migration are complex and, to a large extent, their details still remain unknown. Due to the great depth of subduction the fluid released may be in supercritical condition with, as mentioned earlier, fractal density with strong solvent strength facilitating the hydration and metasomatism of mantle rocks. When the pressure and temperature are reduced to around the critical point, the system goes through a great reduction of gradient of density, accordingly increasing the specific volume which further enlarges porous space and fractures rocks thus in turn facilitating the formation of magma and earthquakes through positive feedback processes.

Self-organized Criticality
The phenomena associated with continuous phase transitions are called critical phenomena, and these are often related to so-called self-organized criticality (SOC). SOC is commonly illustrated conceptually with avalanches resulting from piles of sand which generate a power-law number-size distribution of avalanche magnitudes (Bak et al. 1987). At the criticality point in a SOC phenomenon a small continuous input to the system can cause sudden and discontinuous outputs or avalanches. For example, a fault occurs in broken brittle rock strata when an extra stress is added to change the system at the criticality point. The size and number of faults generated may follow a power law distribution with a small number of large faults and a large number of small faults. SOC is similar to critical point phase transition since both processes involve anomalous state change caused by a minor continuous input pulse at the critical condition point. Numerous studies have also pointed out the effect of the 660-km endothermic phase transition on convection. This could actually generate the periodic occurrence of abrupt changes in convective mode (660-km layered/whole mantle), consecutive with the sudden flushing of oceanic plates previously accumulated above the transition zone (e.g., Le Bars and Davaille 2004).
Many numerical simulations have demonstrated multiple scale and sizeable whole mantle convection, and sublithospheric convection can bring up dense fertile mantle materials from the lower mantle to the upper mantle (Korenaga 2004). Cold downwellings are temporarily stopped by the 660 km endothermic phase change but sink rapidly into the lower mantle (Tackley et al. 1993). The intermittence of layering reflects accumulation and release of negative buoyancy above the endothermic phase boundary (Machetel and Weber 1991;Tackley et al. 1993). The exchange of mass between upper and lower layers can occur in short bursts (Zhong and Gurnis 1994). Although these types of avalanching behaviors are not as easy to test as those of sand piles, one might reasonably assume that these types of processes with SOC nature can generate end products with power law distributions. As a matter of fact, SOC phenomena have been commonly considered to describe extreme geo-events in plate tectonics. Such examples may include but are not limited to earthquakes (Gutenberg and Richter 1944;Turcotte 1997), volcanic eruption durations (Cannavò and Nunnari 2016), plate sizes (Sornette and Pisarenko 2003), slab breakoff (Condie 1998), areal size of magmatism (Pelletier 1999), mineral deposits (Agterberg 1995;Cheng 1999b;Maier and Groves 2011), heat flow over mid-ocean ridges (Cheng 2016), episodic evolution of supercontinents and crustal growth (Cheng 2017b), and energy-probability of earthquakes (Cheng and Sun 2017). Other examples can be found in the book authored by Sornette (2004). The processes involved in response to the preceding extreme events create end products which can be described by frequency-size or frequency-time power law relations. Based on the above reasoning, we may expect lithospheric root detachments and slab breakoffs that occurred during subduction are of difference sizes which follow power-law distributions. Some of these small-sized events may not be noticeable on the surface due to small impact on the global system, but the large detachments and slab breakoffs can cause significant impact on syn-to post-collisional magmatism and metamorphism. The size-frequency distribution of these types of events can be modelled by the following general power-law relation where A represents the size of event and N(>A) the cumulative number of events with size greater than the threshold A. This power-law function involves two constant values: c and b. For example, the well-known Gutenberg-Richter power-law distribution relates the number of large earthquakes to their sizes (Gutenberg and Richter 1944;Turcotte 1997). The exponent, b-value, has been commonly used for predictive purposes. The exponential b-value was found to be internally related to singularity in terms of fractal probability density (Cheng and Sun 2017) with where E(<P) represents the minimum energy released by large earthquakes, with occurrence probability less than P. This equation indicates that the minimum energy released by large earthquakes follows a power-law relation ðβ = 2 3 bÞ for probability of earthquake occurrence with energy greater than E. This model implies that the smaller the probability (P) of a large earthquake, the larger its energy release (E).

Multiplicative Cascade Processes
Multiplicative cascade processes (MCP) are iterative multiplicative processes across multiple scales, which involve positive or negative feedback to generate extreme values that follow multifractal power-law distributions (power-law distributions with multiple exponents) with self-similarities and singularities (Meakin 1987;Scherzter and Lovejoy 1987;Agterberg 2007;Cheng 2014). Examples of MCP are common in the study of geocomplexity such as formation of clouds, severe weather and storms (Scherzter and Lovejoy 1987;Malamud et al. 1996;Turcotte 1997;Veneziano and Furcolo 2002), to just name a few. In terms of mantle convection, the convection processes can be viewed as multiplicative cascade processes that create heterogeneity of the mantle by recycling the materials from upper crust to mantle. On a large scale, Wilson cycle cascade evolution involves the opening and closing of an individual oceanic basin, plate drift, plate subduction and plate collision, involving the recycling of lithosphere material and causing extreme events at the interface of phase transition zones or zones around plate boundaries. Depending on the properties of subduction and other factors, plate subduction may cause slab deformation, erosion and breakoff, deep subduction, and collision of continents. These events are responsible for formation of extreme events such as magmatism and earthquakes. During such processes changes of pressure and temperature as well as water content often provides a positive feedback effect on causes of melting or partial melting of lithosphere and the generation of magma reservoirs and seismicity. In the context of multiplicative cascade processes, the mass and energy distribution resulting from these processes often are proved to have self-similarity and singularity which can be modelled by multifractal distributions (Meakin 1987;Schertzer and Lovejoy 1987;Cheng and Agterberg 2009).
The aforementioned mechanisms (PT, SOC and MCP) can coexist in the evolution of earth dynamics systems which cause cascade effects for anomalous diffusion and strain rate originating earthquakes or magmatism creating flare up formation of magmatic activity or cluster frequency-depth distribution of earthquakes. Based on possible mechanisms (PT, SOC and MCP) corresponding to power-law distributions, the fractal density (power-law density) and the singularity analysis method can be used to characterize the causational relations between extreme events such as magmatic activities and earthquakes and the aforementioned nonlinear mechanisms. In the following section a case study of earthquakes will be used to demonstrate the effect of phase transition on formation and distribution of earthquakes that occur along Pacific plate subduction zones.

Rheology Constitutive Equation
In the study of earth tectonics, rheology is an important concept describing rock properties with respect to flow behavior which can be characterized through the following empirical constitutive equation associating stress and strain rate (e.g., Dimanov et al. 1998).
where ε̇represents the strain rate, σ-the stress, n-the stress exponent; d represents the grain size, m is the grain-size exponent, f H 2 O -the water fugacity, and r-the fugacity exponent, Q-the activation energy, P-the pressure, V-the activation volume, T-the absolute temperature, while R is the molar gas constant, and A-a material constant. The constitutive Eq. (10.21) is often utilized in the literature for describing rheology of ductile crust and since it is so well-known it often is provided without citation and reference. Several authors have investigated this equation by various methods such as by physical experiments (Pharr and Ashby 1983;Dimanov et al. 1998). The parameters involved in the equation can be estimated using a log-linear model except for the last combined term Effects of some of the parameters have been summarized by several authors (e.g., Bürgmann and Dresen 2008). For example, diffusion-controlled deformation is linear in stress with n = 1. Different inverse dependencies on grain size have been predicted for lattice diffusion-and grain boundary diffusion-controlled creep with m = 2 and m = 3, respectively. Creep of fine-grained materials involves grain boundary sliding, which may be controlled by grain boundary diffusion (n = 1) or by dislocation motion (n = 2). For climb-controlled dislocation creep, deformation is commonly assumed to be grainsize insensitive (m = 0) with a stress exponent of n = 3-6 (e.g., Bürgmann and Dresen 2008). Materials for which strain rate is proportional to stress raised to a power n > 1 are referred to as having a power-law rheology, whose effective viscosity (μ = σ ̸ ε̇∝ σ 1 − n ) decreases when stress increases. The significant effects of melt distribution on the rheology of rocks have been reported by many authors (e.g., Dimanov et al. 1998Dimanov et al. , 2000. In general, the strain rate is proportional to the water fugacity. The general bivariate relations between the strain rate and other factors considered in the equation are valid and can be applied to characterize the general associations of factors considered in the system (Wang 2016;Dimanov et al. 2000). However, the equation is valid for normal media that generally do not possess singularity for non-zero values of the factors. It is neither possible to use this equation to describe the singular behaviors of constitutive equation in phase transition nor to directly use it to delineate zones of phase transition. Variable depth-frequency distribution of crustal earthquakes and lithological compositions are often integrated to characterize crust deformation in relation to variations of tectonic styles (Mouthereau and Petit 2003). In the following section my attempt is to derive a proper equation to characterize the rheology in phase transition zones.

Rheology and Phase Transition
In order to explain the phase transition zones in the lithosphere associating the effect of phase transition with origin of seismicity and magmatism, one needs to link the rheology to depth of lithosphere. It has been generally accepted that in the brittle crust, frictional strength increases linearly with depth. Phase transitions separate regions into groups of rocks dominated by quartz, feldspar and olivine, respectively; and these regions are characterized by brittle or plastic properties of lithosphere (e.g., Jackson 2002;Bürgmann and Dresen 2008). It was suggested by Sibson (1974) that brittle strength in the crust can be approximated by the Sibson's formulation in which the coefficients of friction and cohesion for pre-fractured rocks are equal to internal friction and cohesion for intact samples: where σ = σ 1 − σ 3 represents differential stress, z is depth, ρ is average density of the overburden, g is acceleration of gravity, β is a coefficient which depends on the type of faulting, and λ represents the pore fluid ratio. Under hydrostatic pressure, λ is 0.36, and it is 0 and 0.7 for dry and wet conditions, respectively (Mouthereau and Petit 2003). In order to discuss the behavior of rheology around phase transition, let us define depth at the center of the phase transition zone as z 0, which will serve as reference of coordinate for further comparison. Let us also denote a small distance increment (in depth) around the phase transition zone as Δz = abs(z − z 0 ), and the corresponding increment of differential stress around the phase transition zone as Δσ = absfðσ 1 − σ 3 ÞðzÞ − ðσ 1 − σ 3 Þðz 0 Þg. When Δz is very small around the phase transition center z 0 , then we can derive the following approximation assuming changes of depth z, β and λ are neglectable: According to the phase transition property of density and temperature or pressure similar to Eq. (10.17) we can assume the mass density of lithosphere around the phase transition center to be Further assuming that the temperature and depth increments are linearly associated when the depth increment is very small, we obtain Therefore, the derivative of Eq. (10.26) satisfies This result implies that change rate ( Δρ Δz ) of density with depth follows a power-law relation with the increment of depth (Δz). If the exponent b is less than 1, the change rate approaches infinity when Δz → 0, which implies that the change rate of differential stress, according to Eqs. (10.27) and (10.24), can become infinitely large. Assuming the other factors to be negligibly small in Eq. (10.21) when Δz is very small, we obtain If the exponent b is less than 1, then the change of strain rate per increment of depth approaches infinity when Δz → 0. It must be reminded that the derivation of the new Eqs. (10.24-10.28) is based on several assumptions involving first order approximations of factors which may need further theoretical justification (detailed discussion will be published elsewhere). Nevertheless, the results obtained here might be the first power-law model providing possible quantitative description of the singularities of differential stress at the phase transition as indicated in the schematic diagram ( Fig. 10.1).

Frequency-Depth Fractal Density Distribution and Singularity Analysis of Earthquakes
In order to demonstrate the effect of differential stress caused by phase transition on formation and distribution of earthquakes, several datasets of earthquakes with magnitudes three or above were selected for several small regions along the Ring of Fire, the Pacific plate boundaries. Data were downloaded from the USGS website under the section of USGS Earthquake Hazards Program (https://earthquake.usgs. gov/earthquakes/map/). The locations of the 30 small areas selected from Aleutian Islands, Kuril Islands, Mariana, Tonga Trench, Mexico, northern Chile and southern Chile are shown in Fig. 10.2. Several hundreds to thousands of earthquakes are selected in each area. These areas were chosen within a short range from the plate boundaries to ensure they contain enough large earthquakes which occurred along subduction zones with similar properties. The main purpose of the case study here is to validate whether earthquakes that occurred in the subduction zones possess clustering with fractal density; therefore, we choose earthquakes in the depth around the Moho ranging from 30 to 100 km. Considering the issue of depth of shallow earthquakes being set a "normal" depth of 33 km or default depths of 5 or 10 km when depths are poorly constrained by available seismic data, we only analyze the earthquakes with occurring depth ranging 34 to 100 km. The numbers of earthquakes in each dataset were grouped on the basis of 10-km depth frequency bins. A profound peak of frequency distributions can be observed around 33 km in all datasets except for western California. To reduce the effect of the "default peak" at depth 33 km, further analysis of the frequency data will be based on earthquakes with depth from 34 km downward. As an example, the frequency-depth distribution of 1263 earthquakes with magnitude greater or equal to 3 and depths between 34 to 100 km from the Tonga region are shown in Fig. 10.3a with the data grouped in a bin of 10 km (frequency-depth distributions for other datasets are not shown here). This graph shows a profound frequency peak at 34-44 km. By eye examination one can see the frequency around the peak within 60 km (from 34 to 94 km) decaying rapidly from the location of the peak at 34 km downward. To validate the fractal density of frequency clustering distribution, the following local number-depth density of earthquakes around the peak z 0 was constructed where Δz is the window size from z 0 , c and b are two parameters to be estimated using the local singularity analysis method (LSA) with windows of multiple sizes: Δz = 10, 20, …, 60 km. The results are calculated for all 30 datasets. Several selected examples are shown in Fig. 10.3b-h. There is no significant peak at 33 km in the datasets from the areas of western California. The decay curves in Fig. 10.3 are least squares fittings to the data with power-law functions. The results estimated from the six datasets give b = 0.90 (E13), 0.44 (E7), 0.27 (E2), 0.49 (N2), 0.55 (N5), 0.69 (W4) and 0.74 (W11) respectively. Coefficients of determination for the least squares fittings to all six datasets are high with R 2 > 0.98 (student t-value > 14), indicating statistically significant power-law models fitted to the data. The results obtained by local singularity analysis of all 30 datasets (except E1, E3, E8) demonstrate that the frequency-depth distributions for large earthquakes (M ≥ 3) are not uniformly distributed but show clustering which can be modelled by using the local fractal density model of Eq. (10.29). The datasets E1, E3 and E8 show linear decay instead of power-law decay. Moreover, the results (shown as yellow dots in Fig. 10.2) demonstrate that the frequency-depth density distributions of earthquakes from the southwestern boundaries of the Pacific plates depict stronger singularities than those of earthquakes from the southeastern boundaries of Distribution of frequency density of earthquakes with magnitudes equal to or greater than 3 from around Moho at 34 km downward. a Frequency-depth distribution of earthquakes from Tonga region; b-h Distribution of decay of frequency density of earthquakes (#/km) with depths from around peak at 34 km downward; Power-law functions were fitted to the observed data by least squares the Pacific Plates except the earthquakes in conjugation regions of three plate boundaries (e.g., N4, N5, W4, W5, W9-W11, E4, E13, E14) that depict stronger singularity. This finding might be significant for understanding the different mechanisms causing earthquakes between the eastern and western Pacific plate boundaries. As reported in the literature, the western boundaries of the Pacific plates are generally colder and older in comparison with the eastern boundaries (Kong et al. 2016;Okazakl and Hirth 2016). Low slab temperatures resulting from faster subduction cause deeper earthquakes (Wei et al. 2017). Omori et al. (2004) have studied association of the distribution of dehydration events with earthquakes and found non-linear correlation between maximum depth of earthquake and temperature of the slab, with lack of deep earthquakes in young subduction-zones. Their work showed that deeper earthquakes (> 300 km) are mostly located in the selected areas along the western subduction zones of Pacific plates whereas fewer deep earthquakes occurred at the eastern boundaries of Pacific plates. The results of the current research may provide supplementary information about singularity of frequency-depth distribution of shallow earthquakes around Moho in the subductions zones of the Pacific plates. The local singularity analysis may provide a new tool for characterization and distinguishing between earthquakes from a fractal and self-similarity point of view. Further work will extend the analysis to cover more areas and other depths of earthquakes. Other sizes of earthquakes will also be considered.

Discussion and Conclusions
In the first part of the chapter, the purpose of including suggestions about mathematical geosciences or geomathematics as a discipline and introduction to examples of significant contributions of mathematical geoscience scientists to science was to appeal to the public and geoscientists to appreciate the indispensable role that MG can play in the family of geosciences. In the second part of the chapter, the fractal density model was introduced and used for characterizing the power-law rheology of phase transition, and singularity analysis of earthquakes from subduction zones of Pacific plates was demonstrated to be a new and promising nonlinear MG method for modeling extreme and "avalanche" geo-events. Examples of application of singularity analysis not only include earthquakes as introduced in the current chapter but also other types of extreme events such as magmatic flare ups (Cheng 2017a), mid ocean ridge anomalous heat flow (Cheng 2016), flooding caused by tropic storms (Cheng 2008), and mineral deposits as well as ore-caused anomalies in surface media (Cheng 2007). Further comprehensive analysis of earthquakes from other regions and clustering depths will be published in separate papers.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.