This Special Issue of the ANZIAM Journal arose from a conference held at Industrial Research, Gracefield, Lower Hutt, New Zealand, on 29 August 2006, to celebrate the research career of Dr Stephen White. Stephen was an international expert in the development and application of numerical methods to quantify fluid flow in porous media, especially those involving heat and chemical transport. As a mathematical modeller, Stephen White’s work has been published extensively, and applied in many practical areas.

Ignition or thermal explosion in an oxidizing porous body of material can be described by a dimensionless reaction–diffusion equation of the form ∂tu=∇2u+λe−1/u. Here such equations will be formulated in symmetrically shaped bounded regions Ω, effectively reducing the mathematical formulation to that of one dimension. This is critically re-examined from a modern perspective using numerical methods. A computer algorithm is constructed and used to carry out a broad-ranging evaluation of the watershed critical initial temperature conditions for thermal ignition of nonuniform assemblies. It is then shown how the resulting mathematical structure for the ignition threshold curves can be correlated by a hyperbolic conic section with a high degree of accuracy over the full range of positive ambient temperature values. However, this sometimes over-predicts (which is bad) and sometimes under-predicts (which is good) the critical initial condition. The definition of additional dimensionless parameters is found to generate further simplification, leading to a universal correlating form capable of collapsing the entire solution space onto a single line in the plane of the new variables. In addition, this study considers the physically intuitive conjecture that spatial moments of the initial temperature profile ought to possess a direct mathematical link to the critical ignition threshold. As such, the mth-order spatial moment of the critical total energy content integrals is defined, and an empirical result is derived stating that certain orders of this moment should be insensitive to changes in ambient temperature and initial shape profile and may be considered functionally dependent on the dimensionless eigenvalue only, within some quantifiable error band. Spatial moment integrals, based on computed critical threshold conditions, are found to support this conjecture, with the best accuracy obtained for the second-order moments.

The freezing of water to ice is a classic problem in applied mathematics, involving the solution of a diffusion equation with a moving boundary. However, when the water is salty, the transport of salt rejected by ice introduces some interesting twists to the tale. A number of analytic models for the freezing of water are briefly reviewed, ranging from the famous work by Neumann and Stefan in the 1800s, to the mushy zone models coming out of Cambridge and Oxford since the 1980s. The successes and limitations of these models, and remaining modelling issues, are considered in the case of freezing sea-water in the Arctic and Antarctic Oceans. A new, simple model which includes turbulent transport of heat and salt between ice and ocean is introduced and solved analytically, in two different cases—one where turbulence is given by a constant friction velocity, and the other where turbulence is buoyancy-driven and hence depends on ice thickness. Salt is found to play an important role, lowering interface temperatures, increasing oceanic heat flux, and slowing ice growth.

In this paper we look at the performance of trigonometric integrators applied to highly oscillatory differential equations. It is widely known that some of the trigonometric integrators suffer from low-order resonances for particular step sizes. We show here that, in general, trigonometric integrators also suffer from higher-order resonances which can lead to loss of nonlinear stability. We illustrate this with the Fermi–Pasta–Ulam problem, a highly oscillatory Hamiltonian system. We also show that in some cases trigonometric integrators preserve invariant or adiabatic quantities but at the wrong values. We use statistical properties such as time averages to further evaluate the performance of the trigonometric methods and compare the performance with that of the mid-point rule.

Implicit Runge–Kutta methods have a special role in the numerical solution of stiff problems, such as those found by applying the method of lines to the partial differential equations arising in physical modelling. Of particular interest in this paper are the high-order methods based on Gaussian quadrature and the efficiently implementable singly implicit methods.

A three-layer compartmental model of the geological structure in the Taupo Volcanic Zone of New Zealand is developed, based on the assumptions of isostasy (constant geostatic pressure at 25 km depth) and a constant rate of volcanism. The upper layer consists of volcanic infill to a depth of about 2500 m, then a middle layer of greywacke-like material to a depth of about 15 km, and a lower layer of andesitic-like material to a depth of 25 km. Our model assumptions predict that the area of each layer increases at a constant rate; that there is a constant ratio between the rate of energy production from volcanic activity and geothermal convection; and that there is the possibility of an abrupt change from rhyolitic to basaltic volcanism, if the middle layer becomes sufficiently thin. Two models are considered: a rifting and a spreading model. Both models predict the lower layer has an andesitic-like density. The spreading model has difficulty matching heat output with observed extension rates. The rifting model predicts the observed extension rates, but requires very deep circulation of groundwater to be consistent with observed chemical and isotopic properties of geothermal fluids.

Understanding ion transport in conjugated polymers is essential for developing mathematical models of applications of these materials. Previous experimental studies have suggested that cation transport in a conjugated polymer could be either diffusion or drift controlled, with debate over which dominates. In this paper we present a new model of cation transport that explains most of the features seen in a set of recent experiments. This model gives good agreement with measured concentration profiles, except when the profile has penetrated the polymer by more than 60%. The model shows that both diffusion and drift processes can be present. An application of a micro-actuator based on a conjugated polymer is presented to demonstrate that this technology could be used to develop a micro-pump.

This is a review of progress made since [R. McKibbin, “An attempt at modelling hydrothermal eruptions”, Proc. 11th New Zealand Geothermal Workshop 1989 (University of Auckland, 1989), 267–273] began development of a mathematical model for progressive hydrothermal eruptions (as distinct from “blasts”). Early work concentrated on modelling the underground process, while lately some attempts have been made to model the eruption jet and the flight and deposit of ejected material. Conceptually, the model is that of a boiling and expanding two-phase fluid rising through porous rock near the ground surface, with a vertical high-speed jet, dominated volumetrically by the gas phase, ejecting rock particles that are then deposited on the ground near the eruption site. Field observations of eruptions in progress and experimental results from a laboratory-sized model have confirmed the conceptual model. The quantitative models for all parts of the process are based on the fundamental conservation equations of motion and thermodynamics, using a continuum approximation for each of the components.

We present an approximate relation for the effective slip length for flows over mixed-slip surfaces with patterning at the nanoscale, whose minimum slip length is greater than the pattern length scale.

Zinc oxide is known to produce a wide variety of nanostructures that show promise for a number of applications. The use of electrochemical deposition techniques for growing ZnO nanostructures can allow tight control of the morphology of ZnO through the wide range of deposition parameters available. Here we model the growth of the rods under typical electrochemical conditions, using the Nernst–Planck equations in two dimensions to predict the growth rate and morphology of the nanostructures as a function of time. Generally good quantitative and qualitative agreement is found between the model predictions and recent experimental results.

An analytic solution is developed for the one-dimensional dissipational slip gradient equation first described by Gurtin [“On the plasticity of single crystals: free energy, microforces, plastic strain-gradients”, J. Mech. Phys. Solids48 (2000) 989–1036] and then investigated numerically by Anand et al. [“A one-dimensional theory of strain-gradient plasticity: formulation, analysis, numerical results”, J. Mech. Phys. Solids53 (2005) 1798–1826]. However we find that the analytic solution is incompatible with the zero-sliprate boundary condition (“clamped boundary condition”) postulated by these authors, and is in fact excluded by the theory. As a consequence the analytic solution agrees with the numerical results except near the boundary. The equation also admits a series of higher mode solutions where the numerical result corresponds to (a particular case of) the fundamental mode. Anand et al. also established that the one-dimensional dissipational gradients strengthen the material, but this proposition only holds if zero-sliprate boundary conditions can be imposed, which we have shown cannot be done. Hence the possibility remains open that dissipational gradient weakening may also occur.

The current geothermal and volcanic activity in the North Island of New Zealand is explained as a consequence of Pacific and Australian plate interactions over the last 20 million years. The primary hypothesis is that the Kermadec subduction zone has for the last 20 million years or more been retreating in a south-easterly direction at about five centimetres per year. It is surmised that this motion and interaction with another subduction zone almost at right angles to it under the North Island resulted in plate tearing due to the incompatibility of the plate geometry where these subduction zones interacted. The nature and consequences of this plate tearing are partially revealed in published maps of the plate currently under the North Island. If the subducted parts of this plate, as shown in Eiby’s maps, [G. A. Eiby, “The New Zealand sub-crustal rift”, New Zeal. J. Geol. Geophy.7 (1964) 109–133] are straightened, then the plate edge lies on a curve giving a rough picture of their position before being torn and subducted by the Kermadec trench motion. This map of the tear suggests the shape of the edge of a missing plate segment torn from the plate, and implies a rotation of the upper North Island, clockwise approximately 20 degrees, about a point just south of the Thames estuary. A consequence of this plate tearing is that the solid retreating crustal wave generating magma pressure beneath the crest of the solid wave has the potential to inject significant basaltic magma into the crust through the tears. These intrusive magma fluxes have the ability to generate geothermal fields and rhyolitic lavas from crustal melts. This could explain the geothermal activity along the Coromandel peninsula five to seven million years ago, the ignimbrite outcrops about Lake Taupo and the current geothermal and volcanic activity stretching from Taupo to Rotorua.

Mathematical modelling of Wairakei geothermal field is reviewed, both lumped-parameter and distributed-parameter models. In both cases it is found that reliable predictions require five to ten years of history for calibration. With such calibration distributed-parameter models are now used for field management. A prudent model of Wairakei, constructed without such historical data, would underestimate field capacity and provide only general projections of the type of changes in surface activity and subsidence.

There are advantages in viewing orthogonal functions as functions generated by a random variable from a basis set of functions. Let Y be a random variable distributed uniformly on [0,1]. We give two ways of generating the Zernike radial polynomials with parameter l, {Zll+2n(x), n≥0}. The first is using the standard basis {xn,n≥0} and the random variable Y1/(l+1). The second is using the nonstandard basis {xl+2n,n≥0} and the random variable Y1/2. Zernike polynomials are important in the removal of lens aberrations, in characterizing video images with a small number of numbers, and in automatic aircraft identification.