Mathematical Modelling for Earth Sciences Xin-She Yang th e m a t ic a l M o d e l l in g Ma Department of Engineering, University of Cambridge Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 DUNEDIN Published by Dunedin Academic Press ltd Hudson House 8 Albany Street Edinburgh EH1 3QB Scotland www.dunedinacademicpress.co.uk ISBN: 978-1-903765-92-0 © 2008 Xin-She Yang The right of Xin-She Yang to be identiied as the author of this work has been asserted by him in accordance with sections 77 and 78 of the Copyright, Designs and Patents Act 1988 All rights reserved. No part of this publication may be reproduced or transmitted in any form or by any means or stored in any retrieval system of any nature without prior written permission, except for fair dealing under the Copyright, Designs and Patents Act 1988 or in accordance with the terms of a licence issued by the Copyright Licensing Society in respect of photocopying or reprographic reproduction. Full acknowledgment as to author, publisher and source must be given. Application for permission for any other use of copyright material should be made in writing to the publisher. th e While all reasonable attempts have been made to ensure the accuracy of information contained in this publication it is intended for prudent and careful professional and student use and no liability will be accepted by the author or publishers for any loss, damage or injury caused by any errors or omissions herein. This disclaimer does not effect any statutory rights. m a t ic a l M o d e l l in g Ma BrItISH LIBrArY CAtALoguINg IN PuBLICAtIoN DAtA A catalogue record for this book is available from the British Library Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 Printed in the uninted Kingdom by Cpod Contents Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I Mathematical Methods 1 1 Mathematical Modelling 1.1 Introduction . . . . . . . . . . . . . . . . . 1.1.1 Mathematical Modelling . . . . . . 1.1.2 Model Formulation . . . . . . . . . 1.1.3 Parameter Estimation . . . . . . . 1.2 Mathematical Models . . . . . . . . . . . 1.2.1 Differential Equations . . . . . . . 1.2.2 Functional and Integral Equations 1.2.3 Statistical Models . . . . . . . . . 1.3 Numerical Methods . . . . . . . . . . . . . 1.3.1 Numerical Integration . . . . . . . 1.3.2 Numerical Solutions of PDEs . . . 1.4 Topics in This Book . . . . . . . . . . . . E rth 00 for 8) g Ma l in 2 Calculus and Complex Variables 2.1 Calculus . . . . . . . . . . . . . . . . . . 2.1.1 Set Theory . . . . . . . . . . . . 2.1.2 Differentiation and Integration . 2.1.3 Partial Differentiation . . . . . . 2.1.4 Multiple Integrals . . . . . . . . 2.1.5 Jacobian . . . . . . . . . . . . . . 2.2 Complex Variables . . . . . . . . . . . . 2.2.1 Complex Numbers and Functions m a t ic a l M o d e l th e 2.2.2 Analytic Functions . . . . . . . . 2.3 Complex Integrals . . . . . . . . . . . . Xin-She Yang 2.3.1 Cauchy’s Integral Theorem . . . c Dunedin Academic 2.3.2 Residue Theorem . . . . . . . . . a 2 S c ie n c e s ( i iii vii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3 3 5 8 11 11 16 16 17 17 19 20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 23 23 26 33 35 36 38 39 40 41 42 43 th e m a t ic a l M o d e l l in g Ma iiiv Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 CONTENTS ConTenTs 3 Vectors and Matrices 3.1 Vectors . . . . . . . . . . . . . . . 3.1.1 Dot Product and Norm . . 3.1.2 Cross Product . . . . . . . 3.1.3 Differentiation of Vectors . 3.1.4 Line Integral . . . . . . . . 3.1.5 Three Basic Operators . . . 3.1.6 Some Important Theorems 3.2 Matrix Algebra . . . . . . . . . . . 3.2.1 Matrix . . . . . . . . . . . . 3.2.2 Determinant . . . . . . . . 3.2.3 Inverse . . . . . . . . . . . . 3.2.4 Matrix Exponential . . . . 3.2.5 Solution of linear systems . 3.2.6 Gauss-Seidel Iteration . . . 3.3 Tensors . . . . . . . . . . . . . . . 3.3.1 Notations . . . . . . . . . . 3.3.2 Tensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 45 46 47 48 49 49 51 51 51 53 54 54 55 57 58 58 59 4 ODEs and Integral Transforms 4.1 Ordinary Differential Equations . 4.1.1 First-Order ODEs . . . . 4.1.2 Higher-Order ODEs . . . 4.1.3 Linear System . . . . . . 4.1.4 Sturm-Liouville Equation 4.2 Integral Transforms . . . . . . . . 4.2.1 Fourier Series . . . . . . . 4.2.2 Fourier Integral . . . . . . 4.2.3 Fourier Transforms . . . . 4.2.4 Laplace Transforms . . . 4.2.5 Wavelets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 61 62 64 65 66 68 69 73 74 75 77 5 PDEs and Solution Techniques 5.1 Partial Differential Equations . . . . . . . . 5.1.1 First-Order PDEs . . . . . . . . . . 5.1.2 Classification of Second-Order PDEs 5.2 Classic Mathematical Models . . . . . . . . 5.2.1 Laplace’s and Poisson’s Equation . . 5.2.2 Parabolic Equation . . . . . . . . . . 5.2.3 Wave Equation . . . . . . . . . . . . 5.3 Other Mathematical Models . . . . . . . . . 5.3.1 Elastic Wave Equation . . . . . . . . 5.3.2 Reaction-Diffusion Equation . . . . . 5.3.3 Navier-Stokes Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 79 80 80 81 81 82 82 82 83 83 83 . . . . . . . . . . . iiiv CONTENTS ConTenTs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 84 84 87 87 88 89 6 Calculus of Variations 6.1 Euler-Lagrange Equation . . . . . . 6.1.1 Curvature . . . . . . . . . . . 6.1.2 Euler-Lagrange Equation . . 6.2 Variations with Constraints . . . . . 6.3 Variations for Multiple Variables . . 6.4 Integral Equations . . . . . . . . . . 6.4.1 Fredholm Integral Equations 6.4.2 Volterra Integral Equation . . 6.5 Solution of Integral Equations . . . . 6.5.1 Separable Kernels . . . . . . 6.5.2 Volterra Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 91 91 93 99 103 104 104 105 105 105 106 7 Probability 7.1 Randomness and Probability . . . . . . . . 7.2 Conditional Probability . . . . . . . . . . . 7.3 Random Variables and Moments . . . . . . 7.3.1 Random Variables . . . . . . . . . . 7.3.2 Mean and Variance . . . . . . . . . . 7.3.3 Moments and Generating Functions 7.4 Binomial and Poisson Distributions . . . . . 7.4.1 Binomial Distribution . . . . . . . . 7.4.2 Poisson Distribution . . . . . . . . . 7.5 Gaussian Distribution . . . . . . . . . . . . 7.6 Other Distributions . . . . . . . . . . . . . . 7.7 The Central Limit Theorem . . . . . . . . . 7.8 Weibull Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 109 115 116 116 117 118 119 119 120 121 123 124 126 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 131 133 133 133 136 137 137 5.4 5.3.4 Groundwater Flow . . . Solution Techniques . . . . . . 5.4.1 Separation of Variables 5.4.2 Laplace Transform . . . 5.4.3 Fourier Transform . . . 5.4.4 Similarity Solution . . . 5.4.5 Change of Variables . . E rth 00 for 8) g Ma l in 8 Geostatistics 8.1 Sample Mean and Variance . 8.2 Method of Least Squares . . . a 8.2.1 Maximum Likelihood . l c M i ode m at l th e 8.2.2 Linear Regression . . . 8.2.3 Correlation Coefficient Xin-She Yang 8.3 Hypothesis Testing . . . . . . c Dunedin Academic 8.3.1 Confidence Interval . . a 2 S c ie n c e s ( . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv vi CONTENTS ConTenTs 8.4 th e m a t ic a l M o d e l l in g Ma 8.5 Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 8.3.2 Student’s t-distribution . . . . . . . 8.3.3 Student’s t-test . . . . . . . . . . . . Data Interpolation . . . . . . . . . . . . . . 8.4.1 Spline Interpolation . . . . . . . . . 8.4.2 Lagrange Interpolating Polynomials 8.4.3 Bézier Curve . . . . . . . . . . . . . Kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 140 142 142 149 150 151 II Numerical Algorithms 159 9 Numerical Integration 9.1 Root-Finding Algorithms . . . . . 9.1.1 Bisection Method . . . . . . 9.1.2 Newton’s Method . . . . . . 9.1.3 Iteration Method . . . . . . 9.2 Numerical Integration . . . . . . . 9.2.1 Trapezium Rule . . . . . . 9.2.2 Order Notation . . . . . . . 9.2.3 Simpson’s Rule . . . . . . . 9.3 Gaussian Integration . . . . . . . . 9.4 Optimisation . . . . . . . . . . . . 9.4.1 Unconstrained Optimisation 9.4.2 Newton’s Method . . . . . . 9.4.3 Steepest Descent Method . 9.4.4 Constrained Optimisation . . . . . . . . . . . . . . . 161 161 162 164 166 168 168 170 171 173 177 177 178 179 182 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Finite Difference Method 10.1 Integration of ODEs . . . . . . . . . . . 10.1.1 Euler Scheme . . . . . . . . . . . 10.1.2 Leap-Frog Method . . . . . . . . 10.1.3 Runge-Kutta Method . . . . . . 10.2 Hyperbolic Equations . . . . . . . . . . 10.2.1 First-Order Hyperbolic Equation 10.2.2 Second-Order Wave Equation . . 10.3 Parabolic Equation . . . . . . . . . . . . 10.4 Elliptical Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185 185 186 188 188 189 189 190 191 193 11 Finite Volume Method 11.1 Introduction . . . . . . 11.2 Elliptic Equations . . 11.3 Hyperbolic Equations 11.4 Parabolic Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195 195 196 197 198 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v vii CONTENTS ConTenTs 12 Finite Element Method 12.1 Concept of Elements . . . . . . . . . . . . . 12.1.1 Simple Spring Systems . . . . . . . . 12.1.2 Bar Elements . . . . . . . . . . . . . 12.2 Finite Element Formulation . . . . . . . . . 12.2.1 Weak Formulation . . . . . . . . . . 12.2.2 Galerkin Method . . . . . . . . . . . 12.2.3 Shape Functions . . . . . . . . . . . 12.2.4 Estimating Derivatives and Integrals 12.3 Heat Transfer . . . . . . . . . . . . . . . . . 12.3.1 Basic Formulation . . . . . . . . . . 12.3.2 Element-by-Element Assembly . . . 12.3.3 Application of Boundary Conditions 12.4 Transient Problems . . . . . . . . . . . . . . 12.4.1 The Time Dimension . . . . . . . . . 12.4.2 Time-Stepping Schemes . . . . . . . 12.4.3 Travelling Waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III Applications to Earth Sciences 13 Reaction-Diffusion System 13.1 Mineral Reactions . . . . 13.2 Travelling Wave . . . . . . 13.3 Pattern Formation . . . . 13.4 Reaction-Diffusion System . . . . . . . . E rth 00 for 8) g Ma l in 14 Elasticity and Poroelasticity 14.1 Hooke’s Law and Elasticity . 14.2 Shear Stress . . . . . . . . . . 14.3 Equations of Motion . . . . . 14.4 Euler-Bernoulli Beam Theory 14.5 Airy Stress Functions . . . . 14.6 Fracture Mechanics . . . . . . 14.7 Biot’s Theory . . . . . . . . . 14.7.1 Biot’s Poroelasticity . a t ic a l M o d e l ht em 14.7.2 Effective Stress . . . . 14.8 Linear Poroelasticity . . . . . Xin-She Yang 14.8.1 Poroelasticity . . . . . c Dunedin Academic 14.8.2 Equation of Motion . a 2 S c ie n c e s ( . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201 202 202 206 209 209 210 211 215 216 216 218 219 221 221 223 223 225 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 227 227 229 230 231 . . . . . . . . . . . . 235 235 240 241 246 249 252 257 257 259 259 259 262 th e m a t ic a l M o d e l l in g Ma vi viii Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 CONTENTS ConTenTs 15 Flow in Porous Media 15.1 Groundwater Flow . . . . . . . . . . 15.1.1 Porosity . . . . . . . . . . . . 15.1.2 Darcy’s Law . . . . . . . . . 15.1.3 Flow Equations . . . . . . . . 15.2 Pollutant Transport . . . . . . . . . 15.3 Theory of Consolidation . . . . . . . 15.4 Viscous Creep . . . . . . . . . . . . . 15.4.1 Power-Law Creep . . . . . . . 15.4.2 Derivation of creep law . . . 15.5 Hydrofracture . . . . . . . . . . . . . 15.5.1 Hydrofracture . . . . . . . . . 15.5.2 Diagenesis . . . . . . . . . . . 15.5.3 Dyke and Diapir Propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263 263 263 263 265 269 272 277 277 278 283 283 284 285 A Mathematical Formulae A.1 Differentiation and Integration A.1.1 Differentiation . . . . . A.1.2 Integration . . . . . . . A.1.3 Power Series . . . . . . A.1.4 Complex Numbers . . . A.2 Vectors and Matrices . . . . . . A.3 Asymptotic Expansions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 291 291 291 291 292 292 292 293 B Matlab and Octave Programs B.1 Gaussian Quadrature . . . . . B.2 Newton’s Method . . . . . . . B.3 Pattern Formation . . . . . . B.4 Wave Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 295 295 297 299 301 . . . . Bibliography 303 Index 307 Preface E rth 00 for 8) g Ma l in Mathematical modelling and computer simulations are an essential part of the analytical skills for earth scientists. Nowadays, computer simulations based on mathematical models are routinely used to study various geophysical, environmental and geological processes, from geophysics to petroleum engineering, from hydrology to environmental fluid dynamics. The topics in earth sciences are very diverse and the syllabus itself is evolving. From a mathematical modelling point of view, therefore, this is a decision to select topics and limit the number of chapters so that the book remains concise and yet comprehensive enough to include important and interesting topics and popular algorithms. Furthermore, we use a ‘theorem-free’ approach in this book with a balance of formality and practicality. We will increase dozens of worked examples so as to tackle each problem in a step-by-step manner, thus the style will be especially suitable for non-mathematicians, though there are enough topics, such as the calculus of variation and pattern formation, that even mathematicians may find them interesting. This book strives to introduce a wide range of mathematical modelling and numerical techniques, especially for undergraduates and graduates. Topics include vector and matrix analysis, ordinary differential equations, partial differential equations, calculus of variations, integral equations, probability, geostatistics, numerical integration, optimisation, finite difference methods, finite volume methods and finite element methods. Application topics in earth sciences include reaction-diffusion system, elasticity, fracture mechanics, poroelasticity, and flow in porous media. This book can serve as a textbook in mathematical modelling and numerical methods for earth sciences. This book covers many areas of my own research and learning from experts in the field, and it represents my own personal odyssey through the diversity and multidisciplinary exploration. Over these years, I have received valuable help in various ways from my mentors, friends, colleagues, and students. First and foremost, I would like to thank my mentors, tutors and colleagues: A. C. Fowler, C. J. Mcdiarmid and S. Tsou at Oxford University for introducing me to the wonderful world a l c M i ode of applied mathematics; J. M. Lees, C. T. Morley and G. T. Parks at m at l th e Cambridge University for giving me the opportunity to work on the applications of mathematical methods and numerical simulations in Xin-She Yang c Dunedin Academicvarious research projects; and A. C. McIntosh, J. Brindley, K. Seffan and T. Love who have all helped me in various ways. a 2 S c ie n c e s ( ix viii x CONTENTS PReFACe I thank many of my students who have directly and/or indirectly tried some parts of this book and gave their valuable suggestions. Special thanks to Hugo Scott Whittle, Charles Pearson, Ryan Harper, J. H. Tan, Alexander Slinger and Adam Gordon at Cambridge University for their help in proofreading the book. In addition, I am fortunate to have discussed many important topics with many international experts: D. Audet and H. Ockendon at Oxford, J. A. D. Connolly at ETHZ, A. Revil at Colorado, D. L. Turcotte at Cornell, B. Zhou at CSIRO, and E. Holzbecher at WIAS. I would like to thank them for their help. I also would like to thank the staff at Dunedin Academic Press for their kind encouragement, help and professionalism. Special thanks to the publisher’s referees, especially to Oyvind Hammer of the University of Oslo, Norway, for their insightful and detailed comments which have been incorporated in the book. Last but not least, I thank my wife, Helen, and son, Young, for their help and support. While every attempt is made to ensure that the contents of the book are right, it will inevitably contain some errors, which are the responsibility of the author. Feedback and suggestions are welcome. th e m a t ic a l M o d e l l in g Ma Xin-She Yang Cambridge, 2008 Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 Part I th e m a t ic a l M o d e l l in g Ma Mathematical Methods Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 g Ma m a t ic a l M o d e l l in th e Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 Chapter 1 Mathematical Modelling 1.1 1.1.1 Introduction Mathematical Modelling E rth 00 for 8) g Ma l in Mathematical modelling is the process of formulating an abstract model in terms of mathematical language to describe the complex behaviour of a real system. Mathematical models are quantitative models and often expressed in terms of ordinary differential equations and partial differential equations. Mathematical models can also be statistical models, fuzzy logic models and empirical relationships. In fact, any model description using mathematical language can be called a mathematical model. Mathematical modelling is widely used in natural sciences, computing, engineering, meteorology, and of course earth sciences. For example, theoretical physics is essentially all about the modelling of real world processes using several basic principles (such as the conservation of energy, momentum) and a dozen important equations (such as the wave equation, the Schrodinger equation, the Einstein equation). Almost all these equations are partial differential equations. An important feature of mathematical modelling and numerical algorithms concerning earth sciences is its interdisciplinary nature. It involves applied mathematics, computer sciences, earth sciences, and others. Mathematical modelling in combination with scientific computing is an emerging interdisciplinary technology. Many international companies use it to model physical processes, to design new products, a l c M i ode to find solutions to challenging problems, and increase their competim at l th e tiveness in international markets. The basic steps of mathematical modelling can be summarised as Xin-She Yang meta-steps shown in Fig. 1.1. The process typically starts with the c Dunedin Academic analysis of a real world problem so as to extract the fundamental physa 2 S c ie n c e s ( 3 4 Chapter 1. Mathematical Modelling ✘ ✛ Realworld problem ✚ ✙ Physical model (Idealisation) Mathematical model ✛ (PDEs,statistics,etc) ✲ Analysis/Validation (Data, benchmarks) th e m a t ic a l M o d e l l in g Ma Figure 1.1: Mathematical modelling. Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 ical processes by idealisation and various assumptions. Once an idealised physical model is formulated, it can then be translated into the corresponding mathematical model in terms of partial differential equations (PDEs), integral equations, and statistical models. Then, the mathematical model should be investigated in great detail by mathematical analysis (if possible), numerical simulations and other tools so as to make predictions under appropriate conditions. Then, these simulation results and predictions will be validated against the existing models, well-established benchmarks, and experimental data. If the results are satisfactory (which they rarely are at first), then the mathematical model can be accepted. If not, both the physical model and mathematical model will be modified based on the feedback, then the new simulations and prediction will be validated again. After a certain number of iterations of the whole process (often many), a good mathematical model can properly be formulated, which will provide great insight into the real world problem and may also predict the behaviour of the process under study. For any physical problem in earth sciences, for example, there are traditionally two ways to deal with it by either theoretical approaches or field observations and experiments. The theoretical approach in terms of mathematical modelling is an idealisation and simplification of the real problem and the theoretical models often extract the essential or major characteristics of the problem. The mathematical equations obtained even for such over-simplified systems are usually very difficult for mathematical analysis. On the other hand, the field studies and experimental approach is usually expensive if not impractical. Apart from financial and practical limitations, other constraining factors in- 1.1 Introduction 5 clude the inaccessibility of the locations, the range of physical parameters, and time for carrying out various experiments. As computing speed and power have increased dramatically in the last few decades, a practical third way or approach is emerging, which is computational modelling and numerical experimentation based on the mathematical models. It is now widely acknowledged that computational modelling and computer simulations serve as a cost-effective alternative, bridging the gap or complementing the traditional theoretical and experimental approaches to problem solving. Mathematical modelling is essentially an abstract art of formulating the mathematical models from the corresponding real-world problems. The master of this art requires practice and experience, and it is not easy to teach such skills as the style of mathematical modelling largely depends on each person’s own insight, abstraction, type of problems, and experience of dealing with similar problems. Even for the same physical process, different models could be obtained, depending on the emphasis of some part of the process, say, based on your interest in certain quantities in a particular problem, while the same quantities could be viewed as unimportant in other processes and other problems. 1.1.2 Model Formulation E rth 00 for 8) g Ma l in Mathematical modelling often starts with the analysis of the physical process and attempts to make an abstract physical model by idealisation and approximations. From this idealised physical model, we can use the various first principles such as the conservation of mass, momentum, energy and Newton’s law to translate into mathematical equations. Let us look at the example of the diffusion process of sugar in a glass of water. We know that the diffusion of sugar will occur if there is any spatial difference in the sugar concentration. The physical process is complicated and many factors could affect the distribution of sugar concentration in water, including the temperature, stirring, mass of sugar, type of sugar, how you add the sugar, even geometry of the container and others. We can idealise the process by assuming that the temperature is constant (so as to neglect the effect of heat transfer), and that there is no stirring because stirring will affect the effective diffusion coefficient and introduce the advection of water or even vertices in the (turbulent) water flow. We then choose a representative element volume (REV) whose size is very small compared with a l c M i ode the size of the cup so that we can use a single value of concentration to m at l th e represent the sugar content inside this REV (If this REV is too large, there is considerable variation in sugar concentration inside this REV). Xin-She Yang c Dunedin AcademicWe also assume that there is no chemical reaction between sugar and water (otherwise, we are dealing with something else). If you drop a 2 S c ie n c e s ( 6 Chapter 1. Mathematical Modelling ❈ ❈❲ J ✒ Γ Ω REV dS ❍ ❥ ❍ Figure 1.2: Representative element volume (REV). the sugar into the cup from a considerable height, the water inside the glass will splash and thus fluid volume will change, and this becomes a fluid dynamics problem. So we are only interested in the process after the sugar is added and we are not interested in the initial impurity of the water (or only to a certain degree). With these assumptions, the whole process is now idealised as the physical model of the diffusion of sugar in still water at a constant temperature. Now we have to translate this idealised model into a mathematical model, and in the present case, a parabolic partial differential equation or diffusion equation [These terms, if they sound unfamiliar, will be explained in detail in the book]. Let us look at an example. Example 1.1: Let c be the averaged concentration in a representative element volume with a volume dV inside the cup, and let Ω be an arbitrary, imaginary closed volume Ω (much larger than our REV but smaller than the container, see Fig. 1.2). We know that the rate of change of the mass of sugar per unit time inside Ω is t a t ic a l M o d e l h em ∂ ∂t ZZZ cdV, Ω l in g Ma δ1 = Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 where t is time. As the mass is conserved, this change of sugar content in Ω must be supplied in or flow out over the surface Γ = ∂Ω enclosing the region Ω. Let J be the flux through the surface, thus the total mass flux 7 1.1 Introduction through the whole surface Γ is δ2 = ZZ Γ J · dS. Thus the conservation of total mass in Ω requires that δ1 + δ2 = 0, or ∂ ∂t ZZZ Ω cdV + ZZ Γ J · dS = 0. This is essentially the integral form of the mathematical model. Using the Gauss’s theorem (discussed later in this book) ZZ ZZZ J · dS = ∇ · J dV, Γ Ω we can convert the surface integral into a volume integral. We thus have ZZZ ZZZ ∂ ∇ · JdV = 0. cdV + ∂t Ω Ω Since the domain Ω is fixed (independent of t), we can interchange the differentiation and integration in the first term, we now get ZZZ ZZZ ZZZ ∂c ∂c dV + ∇ · JdV = [ + ∇ · J]dV = 0. Ω ∂t Ω Ω ∂t m th e a t ic a l M o d e l where D is the diffusion coefficient which depends on the temperature and the type of materials. The negative sign means the diffusion is opposite to the gradient. Substituting this into the mass conservation, we have l in g Ma Since the enclosed domain Ω is arbitrary, the above equation should be valid for any shape or size of Ω, therefore, the integrand must be zero. We finally have ∂c + ∇ · J = 0. ∂t This is the differential form of the mass conservation. It is a partial differential equation. As we know that diffusion occurs from the higher concentration to lower concentration, the rate of diffusion is proportional to the gradient ∇c of the concentration. The flux J over a unit surface area is given by Fick’s law J = −D∇c, Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s (2 ∂c − ∇ · (D∇c) = 0, ∂t 8 Chapter 1. Mathematical Modelling or ∂c = ∇ · (D∇c). ∂t In the simplified case when D is constant, we have ∂c = D∇2 c, ∂t (1.1) which is the well-known diffusion equation. This equation can be applied to study many phenomena such as heat conduction, pore pressure dissipation, groundwater flow and consolidation if we replace D by the corresponding physical parameters. This will be discussed in greater detail in the related chapters this book. 1.1.3 Parameter Estimation Another important topic in mathematical modelling is the ability to estimate the orders (not the exact numbers) of certain quantities. If we know the order of a quantity and its range of variations, we can choose the right scales to write the mathematical model in the nondimensional form so that the right mathematical methods can be used to tackle the problem. It also helps us to choose more suitable numerical methods to find the solution over the correct scales. The estimations will often give us greater insight into the physical process, resulting in more appropriate mathematical models. For example, if we want to study plate tectonics, what physical scales (forces and thickness of the mantle) would be appropriate? For a given driving force (from thermal convection or pulling in the subduction zone), could we estimate the order of the plate drifting velocity? Of course, the real process is extremely complicated and it is still an ongoing research area. However, let us do some simple (yet not so naive) estimations. t ė = a t ic a l M o d e l h em σ , η l in g Ma Example 1.2: Estimation of plate drifting velocity: we know the drift of the plate is related to the thermal convection, and the deformation is mainly governed by viscous creep (discussed later in this book). The strain rate ė is linked to the driving stress σ by Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 where η is the viscosity of the mantle and can be taken as fixed value η = 1021 Pa s (it depends on temperature). The estimation of η will be discussed in Chapter 15. 9 1.1 Introduction Earth’s surface (Ts ) z ❄ ✻heat flow crust T0 upper mantle Figure 1.3: Estimation of the rate of heat loss on the Earth’s surface. Let L be the typical scale of the mantle, and v be the averaged drifting velocity. Thus, the strain rate can be expressed as ė = v . L Combining this equation with the above creep relationship, we have v= Lσ . η Using the typical values of L ≈ 3000 km ≈ 3 × 106 m, σ ≈ 106 Pa, we have v= 3 × 106 × 106 Lσ ≈ 1.5 × 10−9 m/s ≈ 4.7cm/year. ≈ η 2 × 1021 This value is about right as most plates move in the range of 1 ∼ 10 cm per year. The other interesting thing is that the accurate values of σ and L are not needed as long as their product is about the same as Lσ ≈ 3 × 1012 , the estimation of v will not change much. If we use L ≈ 1000 km ≈ 106 m, then, to produce the same velocity, it requires that σ = 3×106 Pa ≈ 30 atm, or about 30 atmospheric pressures. Surprisingly, the driving stress for such large motion is not huge. The force could be easily supplied by the pulling force (due to density difference) of the subducting slab in the subduction zone. m a t ic a l M o d e l E art 00 for 8) g Ma l in Let us look at another example to estimate the rate of heat loss at the Earth’s surface, and the temperature gradients in the Earth’s crust Xin-She Yang and the atmosphere. We can also show the importance of the sunlight c Dunedin Academic in the heat energy balance of the atmosphere. th e h S c ie n c e s ( 2 10 Chapter 1. Mathematical Modelling Example 1.3: We know that the average temperature at the Earth’s surface is about Ts = 300K, and the thickness of the continental crust varies from d = 35km to 70km. The temperature at the upper lithosphere is estimated about T0 = 900 ∼ 1400K (very crude estimation). Thus the estimated temperature gradient is about T0 − Ts dT = ≈ 9 ∼ 31K/km. dz d The observed values of the temperature gradient around the globe are about 10 to 30 K/km. The estimated thermal conductivity k of rocks is about 1.5 ∼ 4.5 W/m K (ignoring the temperature dependence), we can use k = 3 W/m K as the estimate for the thermal conductivity of the crust. Thus, the rate of heat loss obeys Fourier’s law of conduction q = −k∇T = −k dT ≈ 0.027 ∼ 0.093W/m2, dz which is close to the measured average of about 0.07 W/m2 . For oceanic crust with a thickness of 6 ∼ 7 km, the temperature gradient (and thus rate of heat loss) could be five times higher at the bottom of the ocean, and this heat loss provides a major part of the energy to the ocean so as to keep it from being frozen. If this heat loss goes through the atmosphere, then the energy conservation requires that k dT dT + ka = 0, dz crust dh air where h is the height above the Earth’s surface and ka = 0.020 ∼ 0.025 W/m K is the thermal conductivity of the air (again, ignoring the variations with the temperature). Therefore, the temperature gradient in the air is th e m a t ic a l M o d e l l in g Ma dT k dT =− ≈ −3.6 ∼ −4.5K/km, dh ka dz Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 if we use dT /dz = 30 W/km. The negative sign means the temperature decreases as the height increases. The true temperature gradient in dry air is about 10 K/km in dry air, and 6 ∼ 7K/km in moist air. As the thermal conductivity increases with the humidity, so the gradient decreases with humidity. Alternatively, we know the effective thickness of the atmosphere is about 50 km (if we define it as the thickness of layers containing 99.9% of the air mass). We know there is no definite boundary between the atmosphere and outer space, and the atmosphere can extend up to several 1.2 Mathematical Models 11 hundreds of kilometres. In addition, we can also assume that the temperature in space vacuum is about 4 K and the temperature at the Earth’s surface is 300K, then the temperature gradient in the air is 4 − 300 dT ≈ ≈ −6K/km, dh 50 which is quite close to the true gradient. The higher rate of heat loss (due to higher temperature gradient) means that the heat supplied from the crust is not enough to balance this higher rate. That is where the energy of sunlight comes into play. We can see that estimates of this kind will provide a good insight in the whole process. Of course the choice of typical values is important in order to get a valid estimation. Such choice will depend on the physical process and the scales we are interested in. The right choice will be perfected by expertise and practice. We will give many worked examples like this in this book. 1.2 1.2.1 Mathematical Models Differential Equations E rth 00 for 8) g Ma l in The first step of the mathematical modelling process produces some mathematical equations, often partial differential equations. The next step is to identify the detailed constraints such as the proper boundary conditions and initial conditions so that we can obtain a unique set of solutions. For the sugar diffusion problem discussed earlier, we cannot obtain the exact solution in the actual domain inside the water-filled glass, because we need to know where the sugar cube or grains were initially added. The geometry of the glass also needs to be specified. In fact, this problem needs numerical methods such as finite element methods or finite volume methods. The only possible solution is the long-time behaviour: when t → ∞, we know that the concentration should be uniform c(z, t → ∞) → c∞ (=mass of sugar added/volume of water). You may say that we know this final state even without mathematical equations, so what is the use of the diffusion equation ? The main advantage is that you can calculate the concentration at any time usa l c M i ode ing the mathematical equation with appropriate boundary and initial m at l th e conditions, either by numerical methods in most cases or by mathematical analysis in some very simple cases. Once you know the initial Xin-She Yang c Dunedin Academicand boundary conditions, the whole system history will be determined to a certain degree. The beauty of mathematical models is that many a 2 S c ie n c e s ( 12 Chapter 1. Mathematical Modelling seemingly diverse problems can be reduced to the same mathematical equation. For example, we know that the diffusion problem is governed 2 by the diffusion equation ∂c ∂t = D∇ c. The heat conduction is governed by the heat conduction equation ∂T = κ∇2 T, ∂t κ= K , ρcp (1.2) where T is temperature and κ is the thermal diffusivity. K is thermal conductivity, ρ is the density and cp is the specific heat capacity. Similarly, the dissipation of the pore pressure p in poroelastic media is governed by ∂p = cv ∇2 p, (1.3) ∂t where cv = k/(Sµ) is the consolidation coefficient, k is the permeability of the media, µ is the viscosity of fluid (water), and S is the specific storage coefficient. Mathematically speaking, whether it is concentration, temperature or pore pressure, it is the same dependent variable u. Similarly, it is just a constant κ whether it is the diffusion coefficient D, the thermal diffusivity α or the consolidation coefficient cv . In this sense, the above three equations are identical to the following parabolic partial differential equation ∂u = κ∇2 u. (1.4) ∂t Suppose we want to solve the following problem. For a semi-infinite domain shown in Fig. 1.4, the initial condition (whether temperature or concentration or pore pressure) is u(x, t = 0) = 0. The boundary condition at x = 0 is that u(x = 0, t) = u0 =const at any time t. Now the question what is distribution of u versus x at t? Let us summarise the problem. As this problem is one-dimensional, only the x-axis is involved, and it is time-dependent. So we have ∂2u ∂u = κ 2, ∂t ∂x (1.5) u(x, t = 0) = 0, (1.6) m th e a t ic a l M o d e and the boundary condition l l in g Ma with an initial condition u(x = 0, t) = u0 . (1.7) Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 Let us start to solve this mathematical problem. How should we start and where to start? Well, there are many techniques to solve these 13 1.2 Mathematical Models problems, including the similarity solution technique, Laplace’s transform, Fourier’s transform, separation of variables and others. Similarity variable is an interesting and powerful method because it neatly transforms a partial differential equation (PDE) into an ordinary differential equation (ODE) by introducing a similarity variable ζ, then you can use the standard techniques for solving ODEs to obtain the desired solution. We first define a similar variable ζ= x2 , 4κt (1.8) so that u(x, t) = u(ζ) = f (ζ). Using the chain rules of differentiations ∂ ∂ζ x ∂ ∂ = = , ∂x ∂ζ ∂x 2κt ∂ζ ∂2 x 2 ∂2 1 ∂ 1 ∂ ζ ∂2 = ( + + ) = , ∂x2 2κt ∂ζ 2 2κt ∂ζ κt ∂ζ 2 2κt ∂ζ ∂ ∂ ∂ζ x2 ∂ ζ ∂ = =− =− , ∂t ∂ζ ∂t 4κt2 ∂ζ t ∂ζ (1.9) we can write the PDE (1.5) for u as ζ 1 ′ ζ f ], − f ′ = κ · [ f ′′ + t κt 2κt (1.10) where f ′ = df /dζ. Multiplying both sides by t/ζ, −f ′ = f ′′ (ζ) + 1 ′ f, 2ζ 1 f ′′ = −(1 + ). f′ 2ζ or (1.11) Using (ln f ′ )′ = f ′′ /f ′ and integrating the above equation once, we get ln f ′ = −ζ − 1 ln ζ + C, 2 (1.12) where C is an integration constant. This can be written as (1.13) where K = eC . Integrating it again, we obtain p x a t ic a l M o d e u = f (ζ) = Aerf( ζ) + B = Aerf( √ ) + B, l ht em 4κt (1.14) l in g Ma Ke−ζ f′ = √ , ζ Xin-She Yang where c Dunedin Academic E 00 8) for art h S c ie n c e s (2 2 erf(x) = √ π Z x 2 e−ξ dξ, 0 (1.15) 14 Chapter 1. Mathematical Modelling x ✲ u=u0 u(x, t=0)=0 Figure 1.4: Heat transfer near a dyke through a semi-infinite medium. √ is the error function and ξ is a dummy variable. A = K π and B are constants that can be determined from appropriate boundary conditions. This is the basic solution in the infinite or semi-infinite domain. The solution is generic because we have not used any of the boundary conditions or initial conditions. Example 1.4: For the heat conduction problem near a magma dyke in a semi-infinite domain, we can determine the constants A and B. Let x = 0 be the centre of the rising magma dyke so that its temperature is constant at the temperature u0 of the molten magma, while the temperature at the far field is u = 0 (as we are only interested in the temperature change in this case). The boundary condition at x = 0 requires that Aerf(0) + B = u0 . We know that erf(0) = 0, this means that B = u0 . From the initial condition u(x, t = 0) = 0, we have x ) + u0 = 0. A lim erf( √ t→0 4κt √ Since x/ 4κt → ∞ as t → 0 and erf(∞) = 1, we get A + u0 = 0, or A = −u0 . Thus the solution becomes th e m a t ic a l M o d e l l in x x )] = u0 erfc( √ ), 4κt 4κt where erfc(x) = 1 − erf(x) is the complementary error function. The distribution of u/u0 is shown in Fig. 1.5. g Ma u = u0 [1 − erf( √ Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 From the above solution, we know that the temperature √ variation becomes significant in the region of x = d such that d/ κt ≈ 1 at a 15 1.2 Mathematical Models u/u0 1.0 t = 50 0.5 5 1 0.1 x 0 0 1 2 3 4 5 Figure 1.5: Distribution of u(x, t)/u0 with κ = 0.25. given time t. That is d= √ κt, (1.16) which defines a typical length scale. Alternatively, for a given length scale d of interest, we can estimate the time scale t = τ at which the temperature becomes significant. That is τ= d2 . κ (1.17) This means that it will take four times longer if the size of the hot body d is doubled. Now let us see what it means in our example. We know that the thermal conductivity is K ≈ 3 W/m K for rock, its density is ρ ≈ 2700 Kg/m3 and its specific heat capacity cp ≈ 1000 J/kg K. Thus, the thermal diffusivity of solid rock is κ= K 3 ≈ ≈ 1.1 × 10−6 m2 /s. ρcp 2700 × 1000 (1.18) For d ≈ 1m, the time scale of cooling is τ= 1 d2 ≈ 8.8 × 105 seconds ≈ 10 days. ≈ κ 1.1 × 10−6 (1.19) For a larger hot body d = 100 m, then that time scale is τ = 105 days or 270 years. This estimate of the cooling time scale is based on the el th e assumption that no more heat is supplied. However, in reality, there is usually a vast magma reservoir below to supply hot magma constantly, Xin-She Yang c Dunedin Academicand this means that the cooling time is at the geological time scale over millions of years. a 2 m a t ic a l M o d E 00 8) g Ma l in for rth S c ie n c e s ( 16 1.2.2 Chapter 1. Mathematical Modelling Functional and Integral Equations Though most mathematical models are written as partial different equations, however, sometimes it might be convenient to write them in terms of integral equations, and these integral forms can be discretised to obtained various numerical methods. For example, the Fredholm integral equation can be generally written as Z b u(x) + λ K(x, η)y(η)dη = v(x)y(x), (1.20) a where u(x) and v(x) are known functions of x, and λ is constant. The kernel K(x, η) is also given. The aim is to find the solution y(x). This type of problem can be extremely difficult to solve and analytical solutions exist in only a few very simple cases. We will provide a simple introduction to integral equations later in this book. Sometimes, the problem you are trying to solve does not give a mathematical model in terms of dependent variance such as u which is a function of spatial coordinates (x, y, z) and time t, rather they lead to a functional (or a function of the function u); this kind of problem is often linked to the calculus of variations. For example, finding the shortest path between any given points on the Earth’s surface is a complicated geodesic problem. If we idealise the Earth’s surface as a perfect sphere, then the shortest path joining any two different points is a great circle through both points. How can we prove this is true? Well, the proof is based on the Euler-Lagrange equation of a functional ψ(u) d ∂ψ ∂ψ ), = ( ∂u dx ∂u′ (1.21) where u a function of x, u′ = du/dx, and ψ a function of u(x). The detailed proof will be given later in this book in the chapter dealing with calculus of variations. th e m a t ic a l M o d e l l in g Ma 1.2.3 Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 Statistical Models Both differential equations and integral equations are the mathematical models for continuum systems. Other systems are discrete and different mathematical models are needed, though they could reduce to certain forms of differential equations if some averaging is carried out. On the other hand, many systems have intrinsic randomness, thus the description and proper modelling require statistical models, or to be more specific, geostatistical models in earth sciences. For example, suppose that we carried out some field work and made some observations of a specific quantity, say, density of rocks, over a 17 1.3 Numerical Methods s s s s ❡ A s s s s ❡B s s Figure 1.6: Field observations (marked with •) and interpolation for inaccessible locations (marked with ◦). large area shown in Fig. 1.6. Some locations are physically inaccessible (marked with ◦) and the value at the inaccessible locations can only be estimated. A proper estimation is very important. The question that comes naturally is how to estimate the values at these locations using the observation at other locations? How should we start? As we already have some measured data ρi (i = 1, 2, ..., n), the first sensible thing is to use the sample mean or average of <ρi> as the approximation to the value at the inaccessible locations. If we do this, then any two inaccessible locations will have the same value (because the sample data do not change). This does not help if there are quite a few inaccessible locations. Alternatively, we can use the available observed data to construct a surface by interpolation such as linear or cubic splines. There, different inaccessible locations may have different values, which will provide more information about the region. This is obviously a better estimation than the simple sample mean. Thinking along these lines, can we use the statistical information from the sample data to build a statistical model so that we can get a better estimation? The answer is yes. In geostatistics, this is the well-known Kriging interpolation technique which uses the spatial correlation, or semivariogram, among the observation data to estimate the values at new locations. This will be discussed in detail in the chapter about geostatistics. t l in 1.3.1 g Ma 1.3 a t ic a l M o d e l h em Numerical Methods Numerical Integration Xin-She Yang In the solution (1.14) of problem (1.5), there is a minor problem in the evaluation of the solution u. That is the error function erf(x) because 2 c Dunedin Academic E 00 8) for art h S c ie n c e s ( 18 Chapter 1. Mathematical Modelling it is a special function whose integral cannot be expressed as a simple explicit combination of basic functions, it can only be expressed in terms of a quadrature. In order to get its values, we have to either use approximations or numerical integration. You can see that even with seemingly precise solution of a differential equation, it is quite likely that it may involve some special functions. Let us try to evaluate erf(1). From advanced mathematics, we know its exact value is erf(1) = 0.8427007929..., but how do we calculate it numerically? Example 1.5: In order to estimate erf(1), we first try to use a naive 2 approach by estimating the area under the curve f (x) = √2π e−x in the interval [0, 1] shown in Fig. 1.7. We then divide the interval into 5 equallyspaced thin strips with h = ∆x = xi+1 − xi = 1/5 = 0.2. We have six values of fi = f (xi ) at xi = hi(i = 0, 1, ..., 5), and they are f0 = 1.1284, f1 = 1.084, f2 = 0.9615, f3 = 0.7872, f4 = 0.5950, f5 = 0.4151. Now we can either use the rectangular area under the curve (which underestimates the area) or the area around the curve plus the area under curve (which overestimates the area). Their difference is the tiny area about the curve which could still make some difference. If we use the area under the curve, we have the estimation of the total area as A1 ≈ 0.2(f1 + f2 + f3 + f4 + f5 ) ≈ 0.7686. The other approach gives A2 ≈ 0.2(f0 + f1 + f2 + f3 + f4 ) ≈ 0.91125. Both are about 8% from the true value erf(1) ≈ 0.8247. If we take the average of these two estimates, we get th e A1 + A2 ≈ 0.8399, 2 which is much better, but still 0.3% from the true value. This average method is essentially equivalent to using fi = (fi−1 +fi )/2 to approximate the value of f (x) in each interval. m a t ic a l M o d e l l in g Ma A3 ≈ Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 As you can see from this example, the way you discretise the integrand to estimate the integral numerically can have many variants, subsequently affecting the results significantly. There are much better 19 1.3 Numerical Methods f (x)= √2π e−x −1 0 2 1 Figure 1.7: Naive numerical integration. ways to carry out the numerical integration, notably the Gaussian integration which requires only seven points to get the accuracy of about 9th decimal place or 0.0000001% (see Appendix B). All these techniques will be explained in detail in the part dealing with numerical integration and numerical methods. 1.3.2 Numerical Solutions of PDEs The diffusion equation (1.1) is a relatively simple parabolic equation. If we add a reaction term (source or sink) to this equation, we get the classical reaction-diffusion equation ∂u = D∇2 u + γu(1 − u), ∂t (1.22) th e m a t ic a l M o d e l l in g Ma where u can be concentration and any other quantities. γu(1 − u) is the reaction term and γ is a constant. This seemingly simple partial differential equation is in fact rather complicated for mathematical analysis because the equation is nonlinear due to the term −γu2 . However, numerical technique can be used and it is relatively straightforward to obtain solutions (see the chapter on reaction-diffusion system in this book). This mathematical model can produce intriguing patterns due to its intrinsic instability under appropriate conditions. In the two-dimensional case, we have ∂2u ∂2u ∂u = D( 2 + 2 ) + γu(1 − u). ∂t ∂x ∂y (1.23) Xin-She Yang Using the finite difference method to be introduced in the second half of this book, we can solve this equation on a 2-D domain. Fig. 1.8 2 c Dunedin Academic E 00 8) for art h S c ie n c e s ( 20 Chapter 1. Mathematical Modelling Figure 1.8: Pattern formation of reaction-diffusion equation (1.23) with D = 0.2 and γ = 0.5. shows the stable pattern generated by Eq.(1.23) with D = 0.2 and γ = 0.5. The initial condition is completely random, say, u(x, y, t = 0) =rand(n, n) ∈ [0, 1] where n × n is the size of the grid used in the simulations. The function rand() is a random number generator and all the random numbers are in the range of 0 to 1. We can see that a beautiful and stable pattern forms automatically from an initially random configuration. This pattern formation mechanism has been used to explain many pattern formation phenomena in nature shown in Fig. 1.9, including patterns on zebra skin, tiger skin and sea shell, zebra leaf (green and yellow), and zebra stones. For example, the zebra rocks have reddish-brown and white bands first discovered in Australia. It is believed that the pattern is generated by dissolution and precipitation of mineral bands such as iron oxide as mineral in the fluid percolating through the porous rock. The instability analysis of pattern formation and the numerical method for solving such nonlinear reaction-diffusion system will be discussed in detail later in this book. th e m a t ic a l M o d e l l in g Ma 1.4 Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 Topics in This Book So far, we have presented you with a taster of the diverse topics presented in this book. From a mathematical modelling point of view, the topics in earth sciences are vast, therefore, we have to make a decision to select topics and limit the number of chapters so that the book remains concise and yet comprehensive enough to include important 21 1.4 Topics in This Book (a) (b) (d) (c) (e) Figure 1.9: Pattern formation in nature: (a) zebra skin; (b) tiger skin; (c) sea shell; (c) zebra grass; and (e) zebra stone. E rth 00 for 8) g Ma l in topics and popular numerical algorithms. We use a ‘theorem-free’ approach which is thus informal from the viewpoint of rigorous mathematical analysis. There are two reasons for such an approach: firstly we can focus on presenting the results in a smooth flow, rather than interrupting them by the proof of certain theorems; and secondly we can put more emphasis on developing the analytical skills for building mathematical models and the numerical algorithms for solving mathematical equations. We also provide dozens of worked examples with step-by-step derivations and these examples are very useful in understanding the fundamental principles and to develop basic skills in mathematical modelling. The book is organised into three parts: Part I (mathematical methods), Part II (numerical algorithms), and Part III (applications). In Part I, we present you with the fundamental mathematical methods, including calculus and complex variable (Chapter 2), vector and matrix analysis (Chapter 3), ordinary differential equations and integral transform (Chapter 4), and partial differential equations and classic mathematical models (Chapter 5). We then introduce the calculus of variations and integral equations (Chapter 6). The final two chapters (7 and 8) in Part I are about the probability and geostatistics. In Part II, we first present the root-finding algorithms and numerical integration (Chapter 9), then we move on to study the finite difference and finite volume methods (Chapters 10 and 11), and finite element methods (Chapter 12). In Part III, we discuss the topics as applications in earth sciences. a l We first briefly present the reaction-diffusion system (Chapter 13), then c M i ode m at l th e present in detail the elasticity, fracture mechanics and poroelasticity (Chapter 14). We end this part by discussing flow in porous media Xin-She Yang including groundwater flow and pollutant transport (Chapter 15). c Dunedin Academic There are two appendices at the end of the book. Appendix A a 2 S c ie n c e s ( 22 Chapter 1. Mathematical Modelling is a summary of the mathematical formulae used in this book, and the second appendix provides some programs (Matlab and Octove) so that readers can experiment with them and carry out some numerical simulations. At the end of each chapter, there is a list of references for further reading. References th e m a t ic a l M o d e l l in g Ma Fowler A. C., Mathematical Models in the Applied Sciences, Cambridge University Press, (1997). Gershenfeld N., Nature of Mathematical Modeling, Cambridge University Press, (1998). Kardestruncer H. and Norrie D. H., Finite Element Handbook, McGrawHill, (1987). Kreyszig E., Advanced Engineering Mathematics, 6th Edition, Wiley & Sons, New York, (1988). Murch B. W. and Skinner B. J., Geology Today - Understanding Our Planet, John Wiley & Sons, (2001). Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., Numerical Recipes in C++: The Art of Scientific Computing, 2nd Edition, Cambridge University Press, (2002). Smith G. D., Numerical Solution of Partial Differential Equations, Oxford University Press, (1974). Wang H. F., Theory of Linear Poroelasticity: with applications to geomechanics and hydrogeology, Princeton Univ. Press, (2000). Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 B.2 Newton’s Method 297 if nargin<1, % line 7 help Gauss_quad.m; fstr=‘exp(-x.^2)*2/sqrt(pi)’; a=-1.0; b=1.0; end % line 11 % Converting the input integrand, fstr, to a function f(x) f=inline(fstr,0); % Seven-point integration scheme so zeta_1 to zeta_7 zeta=[-0.9491079123; -0.7415311855; -0.4058451513; 0.0; 0.4058451513; 0.7415311855; 0.9491079123]; % Weighting coefficients w=[0.1294849661; 0.2797053914; 0.3818300505; 0.4179591836; 0.3818300505; 0.2797053914; 0.1294849661]; % Index for the seven points Index=1:7; % line 15 I=(b-a)/2*sum(w(Index).*f((b-a).*(zeta(Index)+1)/2+a)); disp(‘ ’); disp(‘The integral is ’); I % line 16 B.2 Newton’s Method The roots of a function f (x) = 0 can be found using Newton’s iteration method f (xn ) xn+1 = xn − ′ , (B.4) f (xn ) where the initial value xn=0 can be a random guess. So we use x0 =randn where randn is a random number drawn from a normal distribution. In case of multiple roots, it will only produce a single root as the initial guess is random. To get multiple roots, you can call the function many times to obtain all the roots. Newton matlab.m E rth 00 for 8) g Ma l in % Finding roots of f(x)=0 via the Newton’s iteration method % Programmed by X S Yang (Cambridge University) % Usage: Newton(function_str); E.g. Newton(‘x.^5-pi’); % [Notes: Since the initial guess is random, so in case % of multiple roots, only a single root will be given.] function [root]=Newton(fstr) ! line 1 format long; ! line 2 % Default function and values if no input a tica l M ode if nargin<1, l ht em help Newton.m; fstr=‘x.^5-pi’; ! line 5 Xin-She Yang c Dunedin Academicend % Tolerance (to the tenth decimal place) a 2 S c ie n c e s ( 298 Appendix B. Matlab and Octave Programs delta=10^(-10); % Converting the input fstr to a function f(x) f=inline(fstr); % Defining x as the independent variable syms x; % Find the f’(x) from f(x) using diff fprime=inline(char(diff(f(x),x))); % Initial random guess xn=randn; deltax=1; % Iteration until the prescribed accuracy while (deltax>delta) root=xn-f(xn)/fprime(xn); deltax=abs(root-xn); xn=root; end disp(strcat(fstr, ‘ has a root’)); root ! line 10 ! line 15 ! line 17 For example, type in >Newton(‘x.∧5-pi’), it will produces a root root = 1.2572741156, which is accurate to the 10th decimal place. However, as we used some symbolic function in Matlab to obtain f ′ (x), the Octave version is slight different, you have to supply the f ′ (x). For the same example, we now have to type in >Newton octave(’x.∧5-pi’,’5*x.∧4’). th e m a t ic a l M o d e l l in g Ma Newton octave.m Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 % Finding roots of f(x)=0 via the Newton’s iteration method % Programmed by X S Yang (Cambridge University) % Usage: Newton(function,fprime); % E.g. Newton(‘x.^5-3.1415’,‘5*x.^4’); % [Notes: Since the initial guess is random, so in case % of multiple roots, only a single root will be given.] function [root]=Newton(fstr,fprimestr) ! line 1 format long; ! line 2 % Default function and values if no input if nargin<1, help Newton.m; fstr=‘x.^5-3.1415’; fprimestr=‘5*x.^4’; ! line 5 end % Tolerance (to the tenth decimal place) delta=10^(-10); % Converting fstr & fprime to functions f(x) & f’(x) str=strcat(‘(’,fstr); str=strcat(str,‘)/(’); str=strcat(str,fprimestr); str=strcat(str,‘)’); fdivfp=inline(str); ! line 10 % Initial random guess Index Bézier linear, 150 quadratic, 150 1-D, 185, 223 2-D, 230 Airy stress function, 249 algorithms, 186 analytical function, 40 assembly by element, 218 asymptotic error function, 293 Gamma function, 293 Stirling’s formula, 294 autocorrelation, 152 bacteria mobility, 272 bar element, 206 basis function, 148 binomial distribution, 119 Biot’s theory, 257 birthday paradox, 114 bisection method, 162 boundary condition, 222 essential, 219 natural, 219 bulk modulus drained, 258 undrained, 262 E rth 00 for 8) g Ma l in calculus of variations, 16, 91 constraint, 99 curvature, 91 a t ic a l M o d e Dido’s problem, 101 l ht em hanging rope problem, 100 multiple variables, 103 Xin-She Yang pendulum, 99 c Dunedin Academic shortest path, 95 a 2 S c ie n c e s ( central difference, 188 central limit theorem, 124 combination, 113 compaction, 276 complex integral, 41 complex number argument angle, 39 conjugate, 39 Euler’s formula, 40 modulus, 39 complex variables, 39 compressibility, 258 consolidation, 259, 261, 272, 274 consolidation coefficient, 273 constitutive relationship, 276 continuity equation, 265 coordinates cylindrical, 37 polar, 37 spherical, 38 correlation coefficient, 136 covariance, 153 crack propagation, 255 critical point, 26 cross product, 47 cumulative probability function, 121 curl, 49 curvature, 91 Darcy’s law, 263 determinant, 53 diagenesis, 284 diagenetic reaction, 228 differential equation, 11 differentiation, 26 implicit, 29 Leibnitz theorem, 28 partial, 33 307 308 rule, 27 diffusion equation, 11, 266, 267 dilation, 257 divergence, 49 divergence theorem, 292 dot product, 46, 47 DuFort-Frankel scheme, 198 dyke formation, 285 eccentricity, 71 effective stress, 259 elasticity, 235 beam bending, 247 Cauchy-Navier equation, 245 elastostatic, 248 Euler-Bernoulli theory, 246 Hooke’s law, 235 strain tensor, 237 stress tensor, 237 stress-strain relationship, 240 elliptic equation, 193 error function, 18, 168, 177 Euler scheme, 186 Euler-Lagrange equation, 93 exponential distribution, 123 th e m a t ic a l M o d e l l in g Ma finite difference method, 20, 185, 299 finite element method, 201, 219, 301 derivative, 215 Gauss quadrature, 215 finite volume method, 195, 199 fracture energy, 255 fracture mechanics, 251 fracture mode mode I, 252 mode II, 252 mode III, 252 Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 Gauss’s theorem, 51 Gaussian distribution, 121 Gaussian integration, 173, 296 geostatistics, 17, 142, 151 Gibbs free energy, 228 Gneiss, 228 gradient, 49 grain-boundary diffusion, 279 Green’s identity, 51 Green’s theorem, 293 INDEX groundwater flow, 263 Hall-Petch equation, 256 harmonic motion, 99 heat conduction, 81, 192, 198, 223 Hooke’s law, 235, 236 hydraulic conductivity, 266 hydrofracture, 283 hyperbolic equation, 197 first-order, 189 second-order, 190 hypothesis testing, 137 ice ages, 71 inner product, 46 integral multiple, 35 differentiation, 34 Jacobian, 36 integral equation, 104 Fredholm equation, 104 sparable kernel, 105 Volterra equation, 105, 106 integral transform Fourier, 68 Laplace, 75 wavelet, 77 integration, 26, 30 by parts, 31 interpolation Bézier, 142, 150 cubic spline, 144 Lagrange polynomial, 148 linear spline, 142 iteration, 161, 166 iteration method, 57, 194 Gauss-Seidel, 57 J-integral, 256 Jacobian, 36 kriging, 17, 142, 151 Lagrange multiplier, 182 Lagrange polynomial, 148 Lagrangian, 98 Lamé constant, 240 Laplace equation, 81 309 INDEX leap-frog scheme, 188 least-square, 209 linear programming, 184 linear system, 55 load efficiency, 273 log-normal distribution, 124 magma dyke, 14 material derivative, 275 mathematical model, 6, 11 mathematical modelling, 3, 8 matrix, 51 exponential, 54 mean, 117 metaheuristic method, 184 method of least square, 133 Milankovitch cycles, 71 mineral banding, 228 mineral reaction, 227 mineral water, 264 model formulation, 5 moment generating function, 118 Navier-Stokes equation, 83 Newton’s method, 164, 166, 297 Newton-Raphson, 166 normal distribution, 121 nugget effect, 153 null hypothesis, 140 numerical integration, 17, 160, 168 Gauss quadrature, 174, 296 integration point, 175 Simpson’s rule, 172 trapezium rule, 169 numerical method, 19 E rth 00 for 8) g Ma l in obliquity, 71 ODE, 13, 61 optimisation, 177 constrained, 182 hill-climbing, 180 Lagrange multiplier, 182 Newton’s method, 178 a t ic a l M o d e Steepest descent, 179 l ht em unconstrained, 177 Xin-She Yang order notation big O, 170 c Dunedin Academic small o, 171 a 2 S c ie n c e s ( outer product, 47 overpressure, 284 parabolic equation, 191 parameter estimation, 8 pattern formation, 229, 230, 233, 299, 301 instability, 231 PDE, 13, 80, 189 perihelion, 71 permeability, 264 permutation, 113 plate tectonics, 8, 283 Poiseuille flow, 67 Poisson distribution, 119 Poisson’s equation, 216 Poisson’s ratio drained, 259 undrained, 259 pollutant transport, 270 absorption, 271 bacteria, 270 concentration, 271 pore pressure, 257, 266, 273 poroelasticity, 257, 260 porosity, 263, 274 porous media Darcy’s flow, 273 effective pressure, 259 effective stress, 259, 277 groundwater flow, 263 increment of fluid content, 257 isotropic, 257 pollutant transport, 270 pore pressure, 262, 273 pumping test, 267 Terzaghi’s theory, 272 precession, 71 pressure head, 262 probability, 109 axiom, 111 conditional, 113, 115 distribution, 118 event, 109 independent events, 112 median, 118 mode, 118 310 moment, 118 random variable, 110, 116 randomness, 109 sample space, 109 probability density function, 121 th e m a t ic a l M o d e l l in g Ma random variable, 116 continuous, 110 discrete, 110 reaction-diffusion, 227, 231 residue theorem, 42 Riccati equation, 61 root-finding, 161 Runge-Kutta method, 186, 189 Xin-She Yang c Dunedin Academic E 00 8) for art h S c ie n c e s ( 2 scientific computing, 3 sedimentary basins, 284 semivariogram, 17, 152 exponential, 156 Gaussian, 156 linear, 155 spherical, 155 series power, 32 Taylor, 32 set intersect, 24 special, 25 subset, 24 theory, 23 union, 24 shape function, 209 2D, 213 Lagrange polynomial, 212 linear, 211 quadratic, 211 similarity solution, 13 Simpson’s rule, 171 Skempton’s coefficient, 258 specific storage, 258, 265 stability condition, 187, 190 standard normal distribution, 122 stationary point, 26 statistical model, 16 statistics, 131 confidence interval, 137 linear regression, 133 maximum likelihood, 133 INDEX sample mean, 131 sample variance, 131 steady state, 218 Stirling’s series, 171 Stokes’s theorem, 51 strain energy, 256 release rate, 254 stress intensity factor, 251, 252 critical, 253 shape factor, 253 Student’s t-distribution, 138 Student’s t-test, 140 surface energy, 253 tensor, 58 analysis, 59 Cartesian, 59 notations, 58 time-dependent problem, 221 time-stepping, 192 implicit, 187 transient problem, 221 trapezium rule, 169 travelling wave, 229 truss element, 206 uniform distribution, 123 upwind scheme, 190 variance, 117 variogram, 152 vector, 45, 46 vector calculus, 48 Venn diagram, 24, 110 viscosity, 8, 281 viscous creep, 261, 277 void ratio, 263 volume strain, 257 water head, 269 wave equation, 81, 82, 190, 223, 301 weak formulation, 210 weakest link theory, 128 Weibull distribution, 126 Young’s modulus, 237