## Abstract

The paper investigates the probability of failure of two-dimensional and three-dimensional slopes using the random finite-element method (RFEM). In this context, RFEM combines elastoplastic finite-element algorithms with random field theory in a Monte Carlo framework. Full account is taken of local averaging and variance reduction over each element, and an exponentially decaying (Markov) spatial correlation function is incorporated. It is found that two-dimensional probabilistic analysis, which implicitly assumes perfect spatial correlation in the out-of-plane direction, may underestimate the probability of failure of slopes.

## 1. Introduction

A considerable number of studies (e.g. Cavounidis 1987; Duncan 1996; Stark & Eid 1998) have compared the factor of safety from a full three-dimensional slope analysis (FS_{3}) with that obtained from a traditional two-dimensional analysis (FS_{2}) and concluded that in the majority of cases for rather uniform slopes FS_{3}≥FS_{2}. The additional stability observed in three-dimensional slopes is generally attributed to the support offered by the boundary conditions in the out-of-plane direction, which is in contrast to the ‘smooth’ conditions implied in a two-dimensional plane strain analysis. The assumption that two-dimensional analysis leads to conservative factors of safety needs some qualification, however. First, a conservative result will only be obtained if the most ‘pessimistic’ plane within a three-dimensional problem is selected for two-dimensional analysis. Griffiths & Marquez (2008) clearly showed planes within a three-dimensional slope that gave higher (unconservative) two-dimensional factors of safety than that of the full three-dimensional problem. In a slope that contains layering and strength variability in the out-of-plane direction, the choice of a two-dimensional pessimistic section may not be intuitively obvious. Although Hutchinson & Sarma (1985) and Hungr (1987) have both asserted that the factor of safety in three dimension is always greater than that in two dimension, it cannot be ruled out that an unusual combination of soil properties and geometry could lead to a three-dimensional mechanism that is more critical. Bromhead & Martin (2004) argued that some landslide configurations with highly variable cross-sections could lead to failure modes in which the three-dimensional mechanism was the most critical. Other investigators have indicated more critical three-dimensional factors of safety (e.g. Chen & Chameau 1982; Seed *et al*. 1990), although this remains a controversial topic. Most importantly, as the two-dimensional factor of safety is generally considered to be conservative, practitioners are reluctant to invest in the more time-consuming three-dimensional approaches.

Furthermore, it is well known that ‘high’ factors of safety do not necessarily mean low probabilities of failure (e.g. Christian *et al*. 1994; Chowdhury & Xu 1995; Duncan 2000). A key question to be addressed in this paper is: under what circumstances will the probability of failure of a slope predicted by a full three-dimensional analysis be higher than that obtained from an equivalent two-dimensional analysis?

The level of stability of natural and constructed slopes is usually expressed by a *factor of safety*, defined as the ratio of the integral of characteristic shear strength to driving forces (gravitational) over the critical failure surface. In foundations engineering, recent interest and application of load and resistance factor design methods allows the engineer to implicitly account for uncertainties by choosing conservatively high characteristic loads and conservatively low characteristic resistances (e.g. AASHTO 2007). The choice, however, is somewhat arbitrary, and in slope stability analysis, the main uncertainty lies in the characteristic shear strength, which may also depend on groundwater conditions. Slopes with nominally the same factor of safety based on characteristic shear strengths could have significantly different failure probabilities because of the uncertainties and how they are dealt with. Duncan (2000) pointed out that through regulation or tradition, the same value of safety factor is often applied to conditions that involve widely varying degrees of uncertainty. This is not logical. A low safety factor does not necessarily correspond to a high probability of failure and vice versa. The relationship between the factor of safety and probability of failure depends on the uncertainties in loads and resistances.

Finite-element investigations of slope reliability have usually been based on two-dimensional analyses (Paice & Griffiths 1997; Griffiths & Fenton 2000, 2004; Hicks & Samy 2002) with none, to our knowledge, using three-dimensional finite elements. A difficulty with three-dimensional slope analysis is that the methods have tended to be extensions of two-dimensional limit equilibrium methods of slices (to ‘methods of columns’), where soil variability and boundary conditions are hard to account for in a systematic way (e.g. Griffiths *et al*. 2009). Important early work was reported by Vanmarcke (1977), which led other investigators such as Yücemen & Al-Homoud (1990) and Auvinet & Gonzalez (2000) to consider the three-dimensional slope reliability.

A methodology developed by the authors, called the ‘random finite-element method’ (RFEM), is used in this paper for three-dimensional probabilistic analysis (Fenton & Griffiths 2008). The method combines nonlinear finite-element methods with random field generation techniques. This method fully accounts for spatial correlation and averaging and is also a powerful slope stability analysis tool that does not require *a priori* assumptions related to the shape or location of the failure mechanism. In this study, the RFEM is further developed to combine three-dimensional elastoplastic finite elements and three-dimensional random field theory in a Monte Carlo framework to directly assess the influence of the coefficient of variation of soil strength and spatial correlation length on slope reliability. The three-dimensional results are compared with an equivalent two-dimensional probabilistic analysis by RFEM (e.g. Griffiths & Fenton 2004) which assumes plane strain conditions and hence perfect correlation in the out-of-plane direction. It will be shown that under some conditions, three-dimensional slope stability analysis leads to higher probabilities of failure than two-dimensional (see table 1 for the notation used in this paper).

## 2. Deterministic analyses

Before turning to probabilistic analyses, we initially present some deterministic slope stability analysis results involving three-dimensional slopes with uniform (constant) soil properties.

The two-dimensional slope profile shown in figure 1 uses eight-node plane strain finite elements to model a 2h:1v undrained clay slope sitting on rock, with strength parameters *ϕ*_{u}=0 and *C*_{u}=*c*_{u}/(*γ*_{sat}*H*)=0.167. The slope stability analyses use an elastic-perfectly plastic stress–strain with a Tresca failure criterion. The left side of the slope is constrained by vertical rollers and the bottom of the slope is fixed. Using typical slope stability analysis method (stability charts, limit analysis or finite elements), it can be shown that FS≈1.39. A finite-element strength reduction approach (e.g. Griffiths & Lane 1999) gives the nodal displacement vectors at failure shown in figure 2 indicating the general shape and location of the failing soil mass.

A three-dimensional slope profile, modelled using 20-node hexahedral elements in which the cross section of the slope shown in figure 1 is extended by a length *L*/2 in the *z*-direction, is shown in figure 3. The bottom of the mesh (*y*=−*H*) and the side (*z*=0) are fully fixed corresponding to ‘rough’ conditions, while the back (*x*=0) and centre line (*z*=*L*/2) are allowed to move only in a vertical plane (owing to symmetry only half of the total slope length needs to be modelled). The total length *L* of the slope was varied in the range 0.8<*L*/*H*<12, enabling an investigation to be made of the influence of three-dimensionality. Figure 4 shows a typical deformed mesh at failure. A comparison of the factor of safety obtained in the three-dimensional and two-dimensional (plane strain) analyses is given in figure 5. The factor of safety in three dimension was always higher than in two dimension owing to side support, but tended to the plane strain solution for length ratios of the order of *L*/*H*>10. It should be mentioned that Griffiths & Marquez (2007) used a mesh density similar to that shown in figure 3 but demonstrated that while finer meshes always gave slightly lower factors of safety than the coarser mesh, the difference never exceeded 2 per cent.

## 3. Random finite-element method

In the remainder of this study, the (dimensionless) shear strength parameter *C*_{u} is assumed to be a random variable characterized statistically by a lognormal distributions (i.e. the logarithm of the property is normally distributed). The lognormal distribution is one of many possible choices (e.g. Fenton & Griffiths 2008); however, it offers the advantage of simplicity, in that it is arrived by a simple nonlinear transformation of the classical normal (Gaussian) distribution. Lognormal distributions guarantee that the random variable is always positive.

The lognormally distributed undrained shear strength *C*_{u} is characterized by three parameters, the mean, *μ*_{Cu}, the s.d. *σ*_{Cu} and the spatial correlation length . The variability of *C*_{u} can conveniently be expressed by the dimensionless coefficient of variation defined as
3.1
As the actual undrained shear strength field is lognormally distributed, its logarithm yields an ‘underlying’ normally distributed (or Gaussian) field. The mean and s.d. of the underlying normal distribution of are related to the mean and s.d. of *C*_{u} using the standard transformation formulae (e.g. Fenton & Griffiths 2008)
3.2
and
3.3
and their inverse form
3.4
and
3.5
The spatial correlation length is measured with respect to the underlying normal field; hence the spatial correlation length () describes the distance over which the spatially random values of will tend to be significantly correlated spatially. A large value of will imply a smoothly varying field, while a small value will imply a rapidly varying field.

In this work, an exponentially decaying (Markovian) correlation function of the following form is used
3.6
where *ρ*(*τ*) is the correlation coefficient between properties assigned to two points in the random field separated by an absolute distance *τ*.

In the current study, the spatial correlation length has been non-dimensionalized by dividing it by the height of the embankment *H* (figure 3) and will be expressed in the form
3.7
Figure 6*a*,*b* shows typical two-dimensional random fields of undrained strength corresponding to different spatial correlation lengths. Figure 6*a* shows a relatively low spatial correlation length of *Θ*=0.2 (implying a rapidly varying field) and figure 6*b* shows a relatively high spatial correlation length of *Θ*=2 (implying a smoothly varying field). The figures depict the variation of *c*_{u} that has been mapped onto the finite-element mesh and have been scaled in such a way that dark and light regions depict ‘strong’ and ‘weak’ soil, respectively. Black represents the strongest element and white the weakest in the particular simulation. It should be emphasized that both these shear strength distributions come from the same lognormal distribution (same mean and s.d.) and it is only the spatial correlation length that is different.

The input parameters relating to the mean, s.d. and spatial correlation length are assumed to be defined at the ‘point’ level. While statistics at this resolution are obviously impossible to measure in practice, they represent a fundamental baseline of the inherent soil variability, which can be corrected through local averaging to take account of the sample size. In the context of the RFEM approach, the sample size is the volume of each finite element used to discretize the slope, and each element is assigned a constant property. If the point distribution is normal, local averaging results in a reduced variance but the mean is unaffected. In a lognormal distribution, however, both the mean and the s.d. are reduced by local averaging. Following local averaging (Fenton & Vanmarcke 1990), the adjusted statistics (*μ*_{CuA},*σ*_{CuA}) represent the mean and s.d. of the lognormal field that is actually mapped onto the finite-element mesh. Further details can be found in Griffiths & Fenton (2004).

For each simulation of the Monte Carlo process, the random field is initially generated and mapped onto the finite-element mesh. After application of gravity loads, failure is said to have occurred if the algorithm is unable to converge within a user-defined iteration ceiling and tolerance (e.g. Griffiths & Lane 1999). An inability to converge in this context implies that no stress redistribution can be found that is simultaneously able to satisfy both the Tresca failure criterion and global equilibrium with the applied gravitational loads. Each simulation of the Monte Carlo process involves the same mean, s.d. and spatial correlation length of undrained strength; however, the spatial distribution of the property varies from one simulation to the next. Following the Monte Carlo process, *p*_{f} is easily estimated by dividing the number of failures by the total number of simulations. The analysis has the option of including cross-correlation between properties and anisotropic spatial correlation lengths (e.g. the spatial correlation length in a naturally occurring stratum of soil is often higher in the two horizontal directions than in the vertical). For the sake of simplicity in the current study, the spatial correlation length has been assumed to be isotropic (e.g. the same in all three directions). By comparison, two-dimensional plane strain RFEM analyses such as those shown in figure 6 are highly anisotropic because an infinite spatial correlation length (perfect correlation) is implied in the out-of-plane direction.

## 4. Single random variable probabilistic analysis

Before presenting the results of the RFEM analyses, we will first discuss a simplified probabilistic analysis assuming a single random variable (SRV). In this method, the slope is assumed to be uniform (spatially constant properties) with *C*_{u} selected randomly from a lognormal distribution with a given mean and s.d. SRV probabilistic methods actually imply an infinite spatial correlation length .

As there is only one random variable in an undrained slope analysis and , then FS is also lognormally distributed with *v*_{FS}=*v*_{Cu} and the probability of failure (*p*_{f}) is simply equal to the probability that FS will be less than unity. Quantitatively, this equals the area beneath the probability density function of FS corresponding to FS<1. For the slope shown in figure 1 which has *μ*_{FS}=1.39, if we let *v*_{Cu}=*v*_{FS}=0.5, equations (3.2) and (3.3) give the mean and s.d. of the underlying normal distribution of as and , respectively. The probability of failure is therefore given by
4.1
where *Φ* is the cumulative standard normal distribution function.

## 5. Random finite-element method probabilistic analyses

In all the RFEM analyses that follow, the bottom of the mesh (*y*=−*H*) is fully fixed and the back of the mesh (*x*=0) is allowed to move only in a vertical plane. It is noted that unlike the deterministic study shown previously, there is no symmetry in the RFEM analyses owing to the spatial varying soil properties. In these analyses, both rough and smooth boundary conditions have been considered at the ends of the mesh in the out-of-plane direction (*z*=0 and *L*). In the rough cases, the ends are fully fixed, and in the smooth case, they are allowed to move only in a vertical plane. In this study, it was determined that 2000 simulations of the Monte Carlo process for each parametric group were sufficient to give reliable and reproducible estimates of the probability of failure *p*_{f}. It can be noted that neither the rough nor the smooth vertical boundary conditions are particularly realistic. Real three-dimensional slopes tend to have rough sloping sides as might be observed at the abutments of an earth dam. In this paper, however, we have considered only simple boundary condition in order to focus on the influence of three-dimensional failure mechanisms.

Figures 7–9 show typical failed slopes with different (isotropic) correlation lengths given by *Θ*=0.2, 2.0 and 200.0. The grey scale depicts the undrained strength, although it should be emphasized that each figure represents just one simulation sampled from a suite of 2000 Monte Carlo repetitions. It can be seen that the failure zone, when it occurs, typically involves a greater volume of soil when the spatial correlation length is much smaller or much larger than the slope height.

Figure 8 demonstrates an important characteristic in three-dimensional slope analysis called the ‘preferred’ failure mechanism width *W*. This is the width of the failure mechanism in the *z*-direction that the finite-element analysis ‘seeks out’. Over a suite of Monte Carlo simulations, the average preferred failure mechanism width is called *W*_{crit}. It will be shown that this dimension has a significant influence on three-dimensional slope reliability depending on whether the length of the slope *L* is greater than or less than *W*_{crit}.

### (a) Influence of the out-of-plane dimension

For each case shown in table 2, the length ratio (figure 3) is varied in the range 0.2<*L*/*H*<16 to investigate the influence of three-dimensionality, with results presented in figures 10–15. Cases 2–6 are obtained by changing parameters (in bold) from the initial Case 1. All three-dimensional finite elements in the mesh (apart from some on the slope surface) are cubes of side length 2 m. The column marked ‘FS’ gives the factor of safety that would be obtained from a two-dimensional slope stability analysis on an *x*–*y* plane with a uniform strength set equal to *μ*_{Cu}.

In the case of smooth boundary conditions, the *p*_{f} of one slice (*L*/*H*=0.2) in the three-dimensional analysis is equivalent to that given by a two-dimensional RFEM analysis, as the three-dimensional analysis is essentially replicating plane strain. It is also shown in the smooth case that as *L*/*H* is increased, *p*_{f} initially decreases, reaching a minimum before rising to eventually exceed the two-dimensional value.

For given values of *v*_{Cu} and *Θ*, let us define the critical slope length *L*_{crit} and the critical slope length ratio (*L*/*H*)_{crit} as being that value of *L*/*H* for which the slope is safest and its probability of failure *p*_{f} a minimum. It will be shown that this minimum probability of failure in the smooth case occurs when *L*_{crit}≈*W*_{crit}.

If we reduce the slope length ratio below this critical value (*L*<*L*_{crit}), the slope finds it easier to form a global mechanism spanning the entire width of the mesh with smooth end conditions, so the value of *p*_{f} increases, tending eventually to the plane strain value. However, if we increase the slope length ratio above this critical value (*L*>*L*_{crit}), the slope finds it easier to form a local mechanism. As *L*>*W*_{crit}, the mechanism has more opportunities to develop somewhere in the *z*-direction, hence *p*_{f} again increases.

In the rough case, *p*_{f} is close to zero for a narrow slice and increases steadily as *L*/*H* is increased owing to a gradual reduction in the supporting influence of the rough boundaries in the three-dimensional case.

As the length ratio is increased in both the rough and smooth cases, the three-dimensional *p*_{f} eventually exceeds the two-dimensional value, indicating that two-dimensional analysis will always give unconservative results if the slope is long enough. It may also be speculated that *p*_{f}→1 as regardless of boundary conditions.

The critical length ratio (*L*/*H*)_{crit} and the length ratio beyond which the *p*_{f} in a two-dimensional analysis ceases to be conservative (*L*/*H*)_{3>2} are listed in table 3.

When *v*_{Cu} is high, it can be expected that the failure mechanism finds it easier to seek out a localized weak zone, leading to a smaller (*L*/*H*)_{crit}. This is clearly demonstrated if one compares the results of Cases 1 and 2, where increasing the strength variability has reduced not only (*L*/*H*)_{crit} but also the value of the cross-over length ratio (*L*/*H*)_{3>2}.

It can also be observed by comparing Cases 1 and 3 that increasing the spatial correlation length has a similar influence. As the spatial correlation length is reduced, as shown, for example, in figure 7, the soil properties change rapidly from point to point, the locally averaged variance of the undrained strength is greatly reduced and the slope tends to behave more like a deterministic slope with uniform (constant) properties throughout. In this case, the preferred failure mechanism occupies a wider region leading to a higher (*L*/*H*)_{crit}. As the spatial correlation length is increased, however, weak or strong soils are more likely to be bunched together, making local mechanisms easier to develop leading to a lower values of (*L*/*H*)_{crit}.

It is interesting to note from figure 13 (Case 4) that when spatial correlation length and soil variability are both increased, there is no obvious minimum in the *p*_{f} versus *L*/*H* plot, although a critical length ratio of about (*L*/*H*)_{crit}≈2 is indicated at the point where *p*_{f} starts to rise. This case indicates that when the spatial correlation length and soil variability are both high and the slope length ratio is shorter than (*L*/*H*)_{crit}, the preferred mechanism for smooth boundary conditions is essentially global and almost equivalent to the two-dimensional plane strain. As shown in table 3, the cross-over length ratios (*L*/*H*)_{3>2} for both rough and smooth boundary conditions in Case 4 were the lowest of any of the cases considered.

Comparing the results of Cases 1 and 5, reducing the factor of safety (based on the mean) to FS=1.11 reduces both (*L*/*H*)_{crit} and (*L*/*H*)_{3>2}. This is to be expected because local failure is more likely to occur when the soil is weak.

Comparing the results of Cases 5 and 6 with Case 1, it is seen that increasing the slope gradient (from 2*h*:1v to 1h:1v) has a similar influence as reducing the FS. Steeper slopes will have lower (*L*/*H*)_{crit} and (*L*/*H*)_{3>2} than flatter slopes.

In the deterministic factor of safety analyses with rough boundary condition (figure 5), the factor of safety in three dimension was always higher than in two dimension but tended to the plane strain solution for length ratios of the order of *L*/*H*≥10. In the probabilistic analyses, the probability of failure in three dimension with rough boundary conditions was initially lower than in two dimension but increases with increasing length ratio to eventually exceed the two-dimensional value at (*L*/*H*)_{3>2}. Although the rough Cases 1–4 have the same factor of safety (based on the mean) and slope angle as in the deterministic analysis from figure 5, the cross-over length ratio in the probabilistic analyses varied quite widely in the range 3<(*L*/*H*)_{3>2}<18.

### (b) Influence of spatial correlation length

In the following, the coefficient of variation of strength, the factor of safety (based on the mean) and the slope angle have been fixed. The spatial correlation length is varied in the range *Θ*={0.125,0.25,…,8} to investigate the influence of spatial correlation length on the *p*_{f} of different slope length ratios. The parameters are shown in table 4 with results shown in figures 16–19. Also included in these figures are the results obtained by the SRV method and two-dimensional RFEM.

For the very short slope shown in figure 16 (Case 7: *L*/*H*=1) where the slope length ratio is smaller than (*L*/*H*)_{crit}, the three-dimensional *p*_{f} is always lower than the two-dimensional value regardless of boundary conditions. The three-dimensional *p*_{f} with rough boundary conditions is much lower than both the three-dimensional *p*_{f} with smooth boundary and the two-dimensional analysis, confirming that boundary support has a strong influence on *p*_{f} for short slopes.

For the longer slope shown in figure 17 (Case 8: *L*/*H*=6), it can be seen that the two-dimensional analysis underestimates *p*_{f} compared with the three-dimensional analysis for *Θ*>1.5 in the smooth case and *Θ*>3.8 in the rough case. As indicated in figure 7, when the spatial correlation length is small, the weak and strong zones of soil are varying rapidly over short distances and the preferred mechanism tends to be global and occupy a wide zone of soil. As the spatial correlation length is increased in figure 8, the preferred mechanism is attracted to local pockets of weak soil and has less width. If the failure mechanism is of the local type, it has more opportunities to fail at different locations along the slope length direction, leading to higher values of *p*_{f}. As the spatial correlation is further increased, the three-dimensional *p*_{f} eventually increases beyond the two-dimensional *p*_{f} which, for a plane strain analysis, assumes an infinite spatial correlation length in the out-of-plane direction.

Further increase in the length ratio as shown in figure 18 (Case 9: *L*/*H*=12) continues this trend with the cross-over point (*L*/*H*)_{3>2} occurring at still smaller values.

Case 8 (*v*_{Cu}=0.5) shown in figure 17 was reanalysed using a higher soil variability (Case 10: *v*_{Cu}=1.0) with the results presented in figure 19. Comparing these results shows that increased soil variability has caused (*L*/*H*)_{3>2} to be reduced to quite low values for both smooth and rough boundary conditions. It can also be noted that *p*_{f} from both three-dimensional RFEM analyses reaches a maximum at around *Θ*≈1.0 indicating a ‘worst-case’ correlation length. It can be argued, therefore, that a worst-case *Θ* must always exist in a three-dimensional RFEM analyses whenever *p*_{f} exceeds the value corresponding to a slope with uniform (constant) properties throughout. This is because when (implying a uniform slope at each simulation), the preferred failure mechanism must be global as shown in figure 9, and the *p*_{f} from the three-dimensional RFEM analyses must finally return to, or even fall below the SRV value.

It is expected that both the two-dimensional and three-dimensional (smooth) RFEM analyses converge on the SRV solution as *Θ* tends to infinity. For the three-dimensional rough case, however, when , the *p*_{f} tends to a value that depends on the length ratio *L*/*H*. If the length ratio is greater than 10, the rough boundary has little influence on *p*_{f} as indicated previously in the deterministic analysis (figure 5) and *p*_{f} will converge on the SRV solution as . If the slope length ratio *L*/*H* is less than 10, however, the rough boundary will result in *p*_{f} tending to a value below the SRV solution as . To confirm these predictions, Cases 7–9 have been reanalysed with a spatial correlation length set to a very large value and the results summarized in table 5.

## 6. Concluding remarks

The paper has investigated the probability of slope failure using both two-dimensional and three-dimensional RFEM probabilistic analyses. The main conclusion is that by implicitly assuming an infinite spatial correlation length in the out-of-plane direction, two-dimensional (plane strain) probability analysis may underestimate the probability of slope failure. This is counter to the usual assumption made in deterministic slope stability analysis that two-dimensional analysis leads to conservative factors of safety compared with three-dimensional owing to the additional support provided by the boundaries in the out-of-plane direction.

The condition under which two-dimensional ceases to be conservative in probabilistic slope analysis depends on several factors and is problem dependent as shown in table 3. As a general trend, the longer the slope in the out-of-plane direction, the more likely it is that the two-dimensional analysis will underestimate the probability of failure. The paper defined a cross-over length ratio (*L*/*H*)_{3>2} above which a two-dimensional analysis could be unconservative The cross-over length ratio (*L*/*H*)_{3>2} was higher when the boundary conditions in the out-of-plane direction were rough rather than smooth. The lowest value of the cross-over length ratio observed for the more realistic rough cases was (*L*/*H*)_{3>2}≈3, which may be used as a conservative upper limit to the safe use of two-dimensional probabilistic analysis.

It was also observed from the three-dimensional analyses, that for slopes with higher length ratios, a worst-case correlation length leads to a maximum value of the probability of failure. This has implications for design where, in the absence of good quality site-specific data, the worst-case spatial correlation length value should be assumed to ensure conservative probabilistic estimates.

The two-dimensional RFEM program is available from the authors’ website at www.mines.edu/~vgriffit/rfem.

## Acknowledgements

The authors wish to acknowledge the support of NSF grant CMS-0408150 on ‘Advanced probabilistic analysis of stability problems in geotechnical engineering’.

## Footnotes

- Received March 31, 2009.
- Accepted June 30, 2009.

- © 2009 The Royal Society