Banded vegetation is a characteristic feature of semi-arid environments. It occurs on gentle slopes, with alternating stripes of vegetation and bare ground running parallel to the contours. A number of mathematical models have been proposed to investigate the mechanisms underlying these patterns, and how they might be affected by changes in environmental conditions. One of the most widely used models is due to Rietkerk and co-workers, and is based on a water redistribution hypothesis, with the key feedback being that the rate of rainwater infiltration into the soil is an increasing function of plant biomass. Here, for the first time, we present a detailed study of the existence and stability of pattern solutions of the Rietkerk model on slopes, using the software package wavetrain (www.ma.hw.ac.uk/wavetrain). Specifically, we calculate the region of the rainfall–migration speed parameter plane in which patterns exist, and the sub-region in which these patterns are stable as solutions of the model partial differential equations. We then perform a detailed simulation-based study of the way in which patterns evolve when the rainfall parameter is slowly varied. This reveals complex behaviour, with sudden jumps in pattern wavelength, and hysteresis; we show that these jumps occur when the contours of constant pattern wavelength leave the parameter region giving stable patterns. Finally, we extend our results to the case in which a diffusion term for surface water is added to the model equations. The parameter regions for pattern existence and stability are relatively insensitive to small or moderate levels of surface water diffusion, but larger diffusion coefficients significantly change the subdivision into stable and unstable patterns.
In 1950, the geologist William MacFadyen published aerial photographs taken over British Somaliland (now northern Somalia) which provided the first documented evidence of striped patterns of vegetation alternating with bare ground . It is now known that self-organized vegetation patterns are a characteristic feature of semi-arid regions in many parts of the world, particularly Africa [2,3], Australia [4,5] and North America [2,6,7]. On slopes, these patterns consist of stripes running parallel to the contours [8–11]. Many different plant types can be involved in this ‘banded vegetation’, and trees often interact with grasses within the bands . Wavelengths are typically in the range 50–300 m (see table 2 of ).
There are no laboratory replicates of banded vegetation, and fieldwork is difficult because of the remoteness and physical harshness of potential study sites. Therefore, most empirical work has been restricted to remote sensing, most recently using satellite images. This has provided much valuable information; however, it has clear limitations, for instance being unable to make predictions about how patterns might change in response to environmental variations. In view of this, and also as a result of the long space and time scales involved in vegetation pattern formation, mathematical models have emerged as a key research tool. Early models involved either cellular automata [13,14] or large systems of coupled ordinary differential equations . However, partial differential equations are now established as the dominant modelling framework. An early model of this type that has been very influential was due to Klausmeier  and is based on the empirical observation that in semi-arid regions, the infiltration of rain water into the soil is positively correlated with vegetation biomass [8,17–19]. This results in greater water availability per unit of biomass, and thus increased plant growth, at higher vegetation densities—a positive feedback that has the potential to generate spatial patterns . Klausmeier's model  explores this potential in a mathematical framework consisting of coupled reaction–diffusion–advection equations for plant biomass and water density. This model and small extensions of it have been explored in very great detail in both simulation-based research [21–24] and analytical studies [25–32].
Although the Klausmeier model is a valuable tool, it is by construction very simplistic. Perhaps its most pronounced simplification is the use of a single water variable. In reality, water dynamics in semi-arid regions are complex. Rainfall contributes directly to surface water, which must then infiltrate into the soil before becoming accessible to plants; and water uptake within the soil is complicated by spatio-temporal variability in rooting depth [33–35]. Most models building on the Klausmeier equations take some account of this complexity by including separate variables for soil and surface water. An important example of this is the Rietkerk model [36,37], which has been used as the basis for many modelling studies of banded vegetation. The model involves three variables: plant biomass P (gm−2), soil water W (mm) and surface water O (mm), which are functions of space x (m) and time t (days). Note that although the model in [36,37] is formulated in two space dimensions, our work is restricted to one dimension. The equations have the form 1.1a 1.1b 1.1c
Table 1 lists the interpretation of the various parameters. Infiltration of surface water into the soil is taken to be an increasing function of plant biomass. As mentioned above, this variation is observed empirically; it results from higher levels of organic matter in the soil, and the alteration in soil structure caused by the increasing density of roots. The uptake of soil water by plants is assumed to have a Michaelis–Menten-type dependence on soil water, and the plant growth rate is assumed to be proportional to this uptake—this is reasonable for a semi-arid environment where water is the limiting resource. Soil water loss will occur due to both evaporation and drainage (leakage), and for simplicity these processes are assumed to be linear functions of water availability and are combined into a single term. Rainfall is assumed to be constant, but note that Guttal & Jayaprakash  have performed a detailed modelling studying on an adapted version of (1.1) that includes seasonal variation in rainfall. The use of linear diffusion to model soil water flow is deliberately simple; more detailed representations of ground water flow are used in the models of von Hardenberg et al.  and Meron et al.  for vegetation patterning. For plant dispersal, diffusion is again used for reasons of simplicity; a more realistic non-local dispersal term is used in the model of Pueyo et al. . The one-dimensional spatial coordinate x runs in the uphill direction, so that the advection term in (1.1c) represents the downhill flow of surface water. Again some subsequent models have used more detailed terms; in particular, in footnote 18 of , Gilad et al. derive a representation of surface water flow using shallow water theory.
Most previous studies using the Rietkerk model (1.1) have considered pattern formation on flat ground, but in this paper we will focus on banded vegetation patterns on slopes. We will vary the rainfall parameter R, and we fix all other parameters at the values given in Rietkerk et al. , which are listed in table 1. The parameter d merits specific mention because Rietkerk et al.  give a range of values: between 0 and 0.5 d−1. This is because the plant loss term includes any herbivory, the extent of which will vary greatly between sites. We fix d in the middle of this range. For the rainfall parameter R, Rietkerk et al.  again give a range of values: between 0 and 3 mm d−1.
The model (1.1) has two spatially uniform steady states: a ‘desert’ state
and a ‘vegetated’ state (Ps, Ws, Os), where 1.2
For the parameter values given in table 1, these two steady states meet in a transcritical bifurcation at R = 1 mm d−1. For R > 1, the vegetated steady state is stable to homogeneous perturbations, whereas the desert steady state is unstable. For R < 1, the desert steady state is stable, whereas the vegetated steady state is unstable; also Ps < 0 for R < 1 so that this state is not ecologically relevant.
Patterned solutions of (1.1) arise for R > 1 mm d−1 when the vegetated steady state is unstable to spatially inhomogeneous perturbations. This is illustrated in figure 1, which shows large time solutions for R = 1.05, 1.1 and 1.25 mm d−1 on a domain of length 500 m with periodic boundary conditions and with initial conditions consisting of small inhomogeneous perturbations applied to (Ps, Ws, Os). In each case, a periodic spatial pattern develops. The patterns are not stationary: rather they move at a constant speed in the uphill direction. Mathematically, this movement is a consequence of the advection term in (1.1c). In the field, such movement is indeed observed in many cases (see  and table 5 of ) and is due to moisture levels being higher on the uphill edge of vegetation bands than on their downhill edge; this is reflected in lower levels of plant death and higher seedling densities [43,44]. However, some field studies also report stationary banded patterns [2,45,46]. This has been attributed to complicating factors including inhibition of seed germination by long-term changes in soil structure in non-vegetated regions , and preferential dispersal of seeds in the downslope direction, due to transport in run-off [47,48].
2. Existence and stability of pattern solutions
The uphill migration of pattern solutions means that they are periodic travelling waves, with the form , and , where z = x − ct and c is the migration speed. Substituting these solution forms into (1.1) gives a fifth-order system of travelling wave ordinary differential equations. Pattern solutions of (1.1) correspond to limit cycle (periodic) solutions of these equations. Note that our restriction to one space dimension means that the patterns we consider would correspond to stripes running parallel to the contours in a two-dimensional setting. The specification in table 1 fixes all of the parameters in the travelling wave equations except two: the rainfall R and the migration speed c. Our initial objective is to determine the region of the R–c parameter plane in which patterns exist. General theory implies that this region will be bounded by segments of three types of loci: Hopf bifurcations, homoclinic (infinite period) solutions and folds in a limit cycle solution branch. We calculated these loci using wavetrain (www.ma.hw.ac.uk/wavetrain) , which is a software package based on numerical continuation  that is specifically designed for the study of periodic travelling wave solutions of partial differential equations.
Figure 2 illustrates the region of the R–c plane in which patterns exist, determined using wavetrain. It extends to values of the migration speed c that are much too large for ecological realism (up to 14 m d−1) and are also very much larger than those observed in simulations such as figure 1, for which we measured speeds less than 0.1 m d−1 for each of the three values of R. This led us to consider the stability of patterns as solutions of (1.1), which can also be studied using wavetrain. Rademacher et al.  proposed a method for determining the stability of periodic travelling waves by numerical continuation of the spectrum, and this is implemented in wavetrain . Using this approach, we found that patterns are unstable in almost all of the region illustrated in figure 2, with stable patterns only occurring for c less than about 0.2 m d−1. Therefore, we focused attention on that part of the parameter region. We found that there were stable patterns for R between about 0.4 and 1.4 mm d−1, and c between about 0.02 and 0.2 m d−1. At the lower boundary of this region of stable patterns, there is a stability change of Eckhaus (sideband) type, meaning that there is a change in the sign of the curvature of the spectrum at the origin. This stability boundary is a locus of Eckhaus points, which can again be traced using wavetrain . The R–c parameter region giving stable patterns is illustrated in figure 3, in which we also plot several contours of fixed pattern wavelength. This division of the parameter plane into regions giving stable and unstable patterns is of key importance in applications, because only stable patterns will be observed as spatially extensive model solutions, although ‘convectively unstable’ patterns can arise as long-term spatio-temporal transients (see §5).
To consider the parameter region illustrated in figure 3 in more detail, it is helpful to consider the solutions as R varies, with the migration speed c fixed. For values of c above about 0.037 m d−1, and also below about 0.015 m d−1, the pattern solution branch connects the Hopf bifurcation at R ≈ 1.25 mm d−1 to the homoclinic (infinite period) solution at the left-hand edge of the parameter region giving patterns. For smaller values of c, the Hopf bifurcation is subcritical, meaning that the solution branch initially proceeds in the direction of increasing R and then folds. The locus of these folds forms the boundary of the parameter region in which patterns exist: to the right of this locus there are no patterns. Figure 4 shows detail near the fold for c = 0.04 m d−1 and includes plots of the spectra for the two different patterns that coexist at a single value of R. The solutions with larger/smaller period are unstable/stable, respectively. Thus as one moves along the solution branch, there is a change in stability at the fold.
For c between about 0.015 and 0.037 m d−1, the structure is more complicated. There is then a solution branch connecting two homoclinic loci, which folds just beyond the right-most homoclinic locus: the fold and homoclinic loci are very close and cannot be resolved in figure 3. There is also a separate solution branch that emanates from the Hopf bifurcation at R ≈ 1.25 mm d−1. This branch has a snaking form, with a tortuous variation in the period (figure 5). Snaking solution branches have been studied in great detail for localized patterns in the Swift–Hohenberg equation [53,54]; in that context, additional peaks are added to the pattern pulse as one moves up the twists of the snaking branch. By contrast, the snaking branch shown in figure 5 consists entirely of periodic (i.e. non-localized) solutions. The part of the R–c plane in which this snaking behaviour occurs is clarified by plotting the loci of the initial folds along the branch, which form a closed loop in the R–c plane (figure 6). The snaking behaviour occurs for values of c within this loop. Calculation of the spectra suggests that all of the many solutions along these snaking branches are unstable as solutions of (1.1) (figure 7). From these investigations, we conclude that the snaking branches do not have ecological relevance, and those patterns that are stable for c in this range (0.015–0.037 m d−1) lie on the solution branch connecting the two homoclinic loci.
Finally, we comment that for all values of c, there is also a second solution branch that emanates from the part of the Hopf bifurcation locus running through the middle of the pattern region (at R ≈ 1 mm d−1). This branch is very tightly localized in R and involves only patterns with extremely large periods (see figure 5); it is therefore of no practical relevance.
The basic conclusion of this section is that there is a single stable pattern solution in the shaded region of the parameter plane in figure 3. There are additional patterns in some other parts of this region, but they are all unstable. This information is of considerable value in understanding and interpreting the results of model simulations. To illustrate this, we will consider in §3 the way in which the vegetation patterns predicted by (1.1) change as the rainfall parameter is gradually varied. We will show that there are abrupt changes in pattern wavelength which can be explained by the results in figure 3.
3. Pattern evolution for variable rainfall
The model idealization of constant parameter values is of course a simplification. In reality, vegetation dynamics are strongly affected by environmental changes, and as a case study we consider the response to variations in rainfall. We performed a series of long numerical simulations on a domain with periodic boundary conditions. We chose a wavelength compatible with these boundary conditions, i.e. (domain length)/N for some integer N. As an initial condition, we used a pattern of this wavelength (a mode N pattern). To reduce computation time, we chose an initial value of R such that the pattern is relatively close to the stability boundary, on the stable side; our results are not sensitive to the initial choice of R. We calculated the initial pattern using wavetrain. We solved for 20 000 days, which is about 55 years. We then made a small change in R and solved for another 20 000 days. We repeated this process, recording the pattern speed and wavelength immediately before each change in R. We varied the increments in R slightly in order to improve the clarity of the results. Note that we do not add noise or any other external perturbation when we change the value of R. Typical results are shown in figure 8, for two different initial values of R. As R varies, the pattern remains on the period contour until this crosses the stability boundary, when it jumps to a new value. Further changes in R cause the pattern to remain on the new period contour, even when the variation in R is reversed. Simulations starting with patterns of other mode numbers N give similar results.
These results imply that the pattern wavelength depends not only on environmental parameters, but also on their values at previous times. Hysteresis between patterned vegetation and bare ground has been noted in a number of previous modelling papers [9,39,40]. Our results show a different type of hysteresis, between patterns of different wavelength. This type of history-dependence has been noted previously in simulations of the Klausmeier model for banded vegetation [21,55]. Its occurrence in the more realistic Rietkerk model argues strongly that it should also be expected in real ecosystems with banded vegetation if the patterns arise because of water redistribution, because this mechanism is the common assumption underlying both the Klausmeier and Rietkerk models. However, it should be noted that our results are all for finite domains with periodic boundary conditions, and their robustness to changes in boundary conditions have not been investigated. There are very few data against which the prediction of hysteresis can be tested, because it would require measurements of wavelength at a series of time points spanning several decades for a single study site. There are a number of examples of such data for two different time points, which have been collected to assess the extent of uphill migration. Older studies of this type measure pattern location relative to ground benchmarks [56,57] while most recent research is based on satellite images [2,11]. The only dataset that we are aware of involving multiple wavelength estimates over several decades is in . This paper studied banded vegetation at several sites in Niger on six occasions between 1950 and 1995. Rainfall varied significantly over this period, and this was reflected in variations in the band : interband width ratio, but the wavelength of the patterns remained constant. This is consistent with our simulation results and suggests that larger variations in rainfall would be needed to induce shifts in pattern wavelength. Climate change means that such large fluctuations in rainfall and other environmental parameters are increasingly likely, and the high volume and easy availability of satellite images will make any resulting pattern change easy to detect over the coming years.
Our work shows that plots of parameter regions giving patterns are very helpful for understanding results from model simulations, with the boundary between stable and unstable patterns being particularly significant. We do not know how to predict the new wavelength that develops when the current period contour crosses the stability boundary. However, we have found that the new value depends on the increments made in the rainfall parameter R. As an example, the contour of period 50 m crosses the stability boundary at R = 0.705 mm d−1. We started with a (stable) pattern of wavelength 50 m (mode 10) with R = 0.73 mm d−1, again on a domain of length 500 m with periodic boundary conditions. Decreasing R to 0.61 mm d−1 leads to a change in wavelength, with a mode 3 pattern (wavelength 166.7 m) developing between 5000 and 10 000 time units (days) after the change in R. If instead R is reduced to only 0.64 mm d−1, then the solution still resembles a mode 10 pattern after 20 000 time units (days). This is probably due to the very slow growth of the unstable linear modes; a further decrease in R quickly induces a shift to a mode 3 pattern. Conversely, when a larger change in R is applied, to 0.57 mm d−1, vegetation disappears entirely to give the desert steady state, which is stable for R < 1 mm d−1. This is a significant result for real vegetation patterns, suggesting that sufficiently large changes in the environment can eradicate vegetation even at rainfall levels that are high enough to support it.
4. Model extension: lateral spread of surface water
In (1.1), the flow of surface water is assumed to be due entirely to the slope. The surface water variable O is simply the height (mm) above ground of the surface water layer, and if this varies in space then there will be flow even on flat ground, due to pressure differences. In their original model (which assumed flat ground), HilleRisLambers et al.  included a diffusion term for O as a simple representation of this process. When Rietkerk et al.  applied the model to vegetation dynamics on a slope, they removed the diffusion term and replaced it with the advection term that is in (1.1c), on the basis that downhill flow will be dominant. We now consider the effect of including the diffusion term, as well as the advection term, in the equation for surface water. Thus, we add the term DO∂2O/∂x2 to the right-hand side of (1.1c). Again we restrict attention to one space dimension.
With the other parameters fixed at the values given in table 1, we investigated the existence and stability of pattern solutions in the R–c plane as DO is varied (figure 9). Note that the units of DO are m2 d−1. For DO = 1, there is no visible difference from the DO = 0 case (compare figure 9a with figure 3). Increasing DO to 25 has little effect: just a slight shift in the stability boundary. However, for DO above about 35, there is a significant difference. The region in which periodic travelling waves exist is little changed (at least within the part of the R–c plane we are considering) but the stability boundary changes shape completely. Note in particular that the minimum speed for stable patterns decreases significantly. At large DO, the advection term in the O equation is dominated by lateral diffusion, so that our finding of stable patterns that move very slowly is consistent with Rietkerk et al.'s  original observation of stationary patterns on flat ground.
As DO is increased from 50 to 100, the stability change on one part of the stability boundary changes from Eckhaus to Hopf type, meaning that it is associated with eigenvalues away from the origin (see figure 10). This change has major implications for the way in which patterns respond to changes in rainfall. In §3, we showed that as rainfall varies with DO = 0, patterns remain of constant wavelength until the wavelength contour crosses a stability boundary (of Eckhaus type), when there is an abrupt transition to a new wavelength. Corresponding simulations for DO = 100 show that there is no such abrupt transition when the wavelength contour crosses the stability boundary of Hopf type (figure 11). Instead, the pattern remains approximately periodic with the same spatial wavelength, but becomes oscillatory in time. Of course, the pattern is intrinsically oscillatory in time, being a periodic travelling wave, but there are now additional oscillations of higher frequency. Figure 12a,b shows plots of plant density against space and time when R = 1.15 during the simulation used for figure 11. The time course clearly involves a superposition of different oscillations, and this is clarified by calculation of the power spectrum (figure 12c). The dominant temporal period is that associated with the underlying periodic travelling wave, which is equal to the wavelength divided by the wave speed. The power series also shows a second significant period, which corresponds approximately to the most unstable part of the spectrum for the underlying periodic travelling wave.
Among the various mathematical models that have been proposed for vegetation patterning in semi-arid environments, the Rietkerk model [36,37] is probably the most widely used, both as a research tool for ecologists and as the basis for model extensions. Despite this very widespread usage, all work on this model has to our knowledge been simulation-based, except for the numerical bifurcation diagrams presented in Rietkerk et al.'s  original paper. Here, we have performed a systematic study of pattern existence and stability. We have also shown that these results provide a straightforward explanation for the complex, history-dependent changes in pattern that occur in response to variations in rainfall. Finally, we considered the way in which our results are affected by the inclusion of a diffusion term in the equation for surface water. This showed that abrupt changes in wavelength only occur when the wavelength contour crosses stability boundaries of Eckhaus type. When the surface water diffusion coefficient is sufficiently large, one of the stability boundaries is of Hopf type, and a crossing of this boundary results in a gradual change in pattern form and the onset of more complicated temporal oscillations, rather than an abrupt transition to a different pattern. All of our work is in one space dimension, so that the patterns that we consider correspond to stripes running parallel to the contours. Extension of our work to two space dimensions is an important but challenging goal for future research.
A key aspect of our work is the classification of patterns into those that are stable/unstable as solutions of (1.1). It is important to emphasize that unstable patterns are not necessarily irrelevant for real instances of banded vegetation. Unstable solutions of a partial differential equation subdivide according to whether the instability is ‘convective’ or ‘absolute’ [59–61]. In the former case, the solutions can occur as persistent spatio-temporal transients [62,63]. Unfortunately, the numerical procedures currently available to distinguish convective and absolute instabilities apply only to very special types of partial differential equation that do not include (1.1). However, we can say definitely that both types of instability do occur in the Rietkerk model. General theory implies that patterns sufficiently close to the stability boundary (and on the unstable side) will be convectively unstable . Moreover, if the spectrum of the pattern contains an isola to the right of its unbounded part and in the right-hand half of the complex plane, then it is known that the solution is absolutely unstable; this was proved by Rademacher  and his paper contains a precise statement of the result. Thus, for example, ‘wave 1’ in figure 4 is absolutely unstable, as are all six solutions on the snaking branch shown in figure 7. Numerical methodology for determining absolute stability is an active research area, and if suitable methods become available then a systematic subdivision of the parameter region giving unstable waves would be a valuable study.
Vegetation patterns in semi-arid environments have been studied for more than 60 years. Within the last 10–15 years, it has become clear that they are just one example of self-organized patterns at the landscape scale. A comprehensive survey of such patterns is given in Rietkerk & van de Koppel , and we discuss here a few examples. Many savannah ecosystems comprise localized patches of trees in a grassy background . These patterns arise primarily from interactions between rainfall and fire; the former is often erratic and its impact depends on both soil structure and local topography , while the frequency of fires has a strong negative correlation with tree cover . The result is a complex and dynamically evolving vegetation structure . In peatlands, one commonly observes patterns of ridges and hollows, with the latter sometimes containing water pools [70–72]. The ridges contain a layer of aerobic peat (the ‘acrotelm’) that contains peat-forming plants and mosses, and which is thin or absent in the hollows. There is a positive feedback between acrotelm thickness and the rate of peat formation , and the observed patterns are thought to arise from the combination of this feedback, a scale-dependent accumulation of nutrients in the ridges  and water-ponding due to lower hydraulic conductivity in the ridges compared with the hollows .
In subalpine forests, strong and consistent winds can cause trees to self-organize into linear patterns . The key mechanism here is that one tree shelters another on its downwind side. The resulting linear patterns can run either parallel or perpendicular to the prevailing wind direction. The former pattern type is known as ‘hedges’ and is widespread, having been documented in North America, Japan and New Zealand . Lines of trees running perpendicular to the wind direction and separated by unforested gaps are known as ‘ribbon forests’ ; they appear to be restricted to the Rocky Mountains. Snow drifting along the lines of trees is thought to be an important factor in the persistence of these patterns [79,80]: the deep parts of snow drifts prevent seedling establishment, whereas shallower parts are beneficial since they provide shelter during the winter. Geomorphology is also thought to play an important role in the formation of some ribbon forests, with trees lying along geological ridges . ‘Wave regenerated forests’ (a.k.a. ‘Shimagare’) are a different type of linear tree pattern, also running perpendicular to the prevailing wind direction [82,83]. Here, there are alternating bands of live and dead trees; in the former, tree height and age increases in the upwind direction up to an abrupt interface, while extensive regeneration occurs among the dead and dying trees. The bands gradually move in the windward direction, on the timescale of the tree generation time. Patterns of this type occur in Japan, USA and Argentina . As a final example, we mention mussel beds in the Wadden Sea, a large intertidal region bordering the Netherlands, Germany and Denmark. These are self-organized into stripes perpendicular to the tidal flow, with a wavelength of 6–10 m . Two different mechanisms have been proposed for this type of pattern formation. Van de Koppel et al.  argue that the binding of mussels to one another, via byssal threads, could cause reduced losses from predation and wave dislodgement at higher mussel densities. Their model has been studied in great detail by Wang et al. . Van Leeuwen et al.  suggest that the greater deposition of sediment under large clumps of mussels could provide an alternative explanation: the increased elevation gives mussels greater access to algal food in higher water layers. Liu et al.  present a detailed model demonstrating the feasibility of this mechanism, and comparing its implications with the reduced losses hypothesis.
Some of these landscape-scale patterns are more amenable to empirical study than banded vegetation in semi-arid landscapes. For example, van de Koppel et al.  successfully achieved self-organized patterning of mussel beds in laboratory tanks, although the patterns were labyrinthine rather than banded. But in all cases, the predictive ability of experiments is severely limited, and mathematical modelling has a vital role to play in understanding how the patterns will be affected by changing environmental conditions. Our work suggests that detailed mathematical investigations of pattern formation in such models will provide a key framework for the understanding of simulation-based studies.
A.S.D. was supported by the Centre for Analysis and Nonlinear PDEs funded by the UK EPSRC grant no. EP/E03635X and the Scottish Funding Council.
- Received May 3, 2014.
- Accepted July 25, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.