The Proceedings of the Eighth International Conference on Creationism (2018)

explicit method T i,j is thus updated as follows: where μ ≡ ᴋ ∆ t/( ∆ z) 2 , a parameter which serves as a dimensionless timestep; ∆t is the actual timestep and ∆z the mesh spacing. Equation (11) is based on a locally quadratic approximation to the temperature field. Thus the (upward) surface heat flux, where j =0, is given by Stability of this scheme requires that μ ≤ 0.5 (Lapidus and Pinder 1982, section 4.4). For the scheme to produce a solution corresponding to the solution of the original heat diffusion equation, the key property of convergence , the original PDE (including boundary and initial conditions) must be well-posed, the difference equation (11) must be consistent with it, and the scheme must be stable. This further implies the necessity for mesh convergence , i.e. as the mesh is refined in both t and z directions, the solution must tend to a final form. In cases involving a heat sink characterized by a function h ( z , t ), the fundamental heat diffusion equation (1) becomes where in our spreadsheets h ( z , t ) is given in units of Wm -3 (watts per cubic metre), ρ is density in kgm -3 and C ρ is specific heat capacity in Jkg -1 K -1 . In practice, the final form of heat sink is made proportional to the difference between the temperature at a point within the plate and the temperature there with the system in thermal equilibrium. The difference equation (11) thus becomes where τ hs is our choice of heat sink lifetime. The inclusion of a heat sink in the difference equation means that the analysis underlying the stability criterion μ ≤ 0.5 given by Lapidus and Pinder (1982) no longer applies. However where we have included a heat sink the timestep is very small (μ < 10 -5 ) and we can use observed mesh convergence to strengthen confidence that our calculations in these cases are essentially correct and stable in that no result depends critically on the timestep and mesh size. Equations (11) and (14) are readily implemented in our spreadsheets, together with the fixed temperature boundary conditions at the top and bottom of the plate and the initial condition simply as the first column of T i,j values with i =0 and j =0, 1, . . . n . (i.e. T i=0,j=0 = T 0 , T i=0,j=1,2, . . .n = T 1 ). The total heat loss is calculated using equation (5), where T m is the average of all the mid-interval temperatures through the lithospheric column, and then the shrinkage from equation (2). However calculation of the surface heat flux is complicated by the existence of a thermal boundary layer near the surface and associated mesh resolution issues. These are considered in the following subsection. 3. The thermal boundary layer and mesh resolution In mathematical terms a boundary layer is a narrow region where the solution of a differential equation changes rapidly. It occurs when the highest-order term in the equation is multiplied by a small parameter ε << 1 and must have the property that the layer thickness δ → 0 as ε → 0 (Bender and Orszag 1978). Here in physical terms we are dealing with the downward propagation of a broadening cooling front from the surface of initially hot material subject to the heat diffusion equation (1). At relatively early times, when the cooling front has only propagated a short distance downwards such that it has not yet felt the influence of the bottom boundary, the solution of the half space model is a good approximation to the cooling process. In this case the temperature field depends only on the combination and after time t the cooling front will have propagated a distance of order This propagation distance defines a thermal boundary layer whose thickness may be formalized as such that equation (3) for the surface heat flux in the half space model may be written in the form If the timestep Δ t is much smaller than the stability limit, i.e., μ << 0.5, as is the case in our short time scale calculations with κ retaining its natural value (≈8×10 -7 m 2 /s), there will be a thin thermal boundary layer, defined by δ th <<Δ z for some time. In these cases the surface heat flux, shrinkage and temperatures through the column of lithosphere are best calculated from the half space model because the spatial discretization is too coarse to resolve the near-surface temperature variation. However this is only possible in the absence of a heat sink since the half space model solution is no longer valid if there is a heat sink in operation. In the long time scale calculation, with Δ z = 950 m and Δ t = 10,000 years, the boundary layer thickness after a single timestep is 893 m, which is comparable with Δ z . However close to the surface the heat flux and temperature variation are not well approximated locally by a quadratic form; equations (11) and (12) still give inaccurate results. In these circumstances the results for the half space model are used instead until after a few hundred timesteps, when the mesh resolution has become adequate and the now very broad cooling front is just about to interact noticeably with the bottom boundary. In our model this switch-over point, although somewhat arbitrary, is taken to be at 8.0 Ma. In the above cases where mesh resolution and the curvature of the temperature profile render the discretized calculation inaccurate, special care has to be taken in the calculation of total heat loss and consequently of shrinkage. In principle the total heat loss is calculated from the sum of heat losses for individual z intervals, which depend on the drop in average temperature for each. The simplest way of calculating this sum is by the trapezoidal rule , which treats the temperature variation between mesh points as linear. For cases with a significantly curved temperature profile we employ Simpson’s rule instead; this is based on a piecewise quadratic approximation to the temperature profile, and for the most extreme points (i.e. adjacent to the plate surface in the first few timesteps) gives an order of magnitude improvement in the accuracy of the result. It is not needed beyond the switch-over Worraker and Ward ◀ Ocean floor cooling ▶ 2018 ICC 676