## Abstract

We derive an analytical expression for entropy production in earthquake populations based on Dewar’s formulation, including flux (tectonic forcing) and source (earthquake population) terms, and apply it to the Olami–Feder–Christensen numerical model for earthquake dynamics. Assuming the commonly observed power-law rheology between driving stress and remote strain rate, we test the hypothesis that maximum entropy production (MEP) is a thermodynamic driver for self-organized ‘criticality’ (SOC) in the model. MEP occurs when the global elastic strain is near-critical, with small relative fluctuations in macroscopic strain energy expressed by a low seismic efficiency, and broad-bandwidth power-law scaling of frequency and rupture area. These phenomena, all as observed in natural earthquake populations, are hallmarks of the broad conceptual definition of SOC (which has, to date, often included self-organizing systems in a near but strictly subcritical state). In the MEP state, the strain field retains some memory of past events, expressed as coherent ‘domains’, implying a degree of predictability, albeit strongly limited in practice by the proximity to criticality and our inability to map the natural stress field at an equivalent resolution to the numerical model.

## 1. Introduction

The maximum entropy principle is useful in determining the underlying distribution of energy states at equilibrium in a variety of physical systems with known macroscopic constraints (Jaynes 1957). In general, this results in a probability density distribution of the Boltzmann form *p*(*E*)∼*g*(*E*)e^{−E/kT}, where *g* is the degeneracy of states, *k* is Boltzmann’s constant and *T* is the temperature. In recent years, it has been suggested that the maximum entropy principle may also hold for some far-from-equilibrium, steady-state, dissipative systems. For example, the twin constraints of finite tectonic energy flux and the fact that source rupture area is confined to a near-two-dimensional fractal set of potential fault rupture surfaces imply a gamma frequency distribution of radiated energy, *E*_{S}, of the form , where *B* is a power-law exponent and *θ* is a characteristic energy (Main & Burton 1984; Main 1996). For self-similar ruptures, the radiated energy is related to the rupture area by *E*_{S}∼*A*^{3/2}, so *p*(*A*)∼*A*^{−b−1}*e*^{−(A/θA)3/2}, where *b*=(3/2)*B* is equal to the Gutenberg–Richter exponent *b* in the frequency (*F*) magnitude (*m*) relation, (Kanamori & Anderson 1975; Turcotte 1997). Additional constraints may be appropriate, for example, angular momentum or enstrophy are commonly used in fluid systems, and may be appropriate in modelling solid-state flow in the whole Earth system, but these are negligible in terms of the brittle field plate tectonic forcing we investigate here.

Real data on regional earthquake populations (Kagan 1997) do commonly show a gamma distribution of scalar seismic moment (the product of the rigidity, rupture area and average slip), a source parameter proportional to the radiated energy. However, it can be hard to pin down *θ*, which depends on the distribution of extreme events that may not have been recorded in current time-limited datasets (Main *et al.* 2008), far less converged to a stable central limit (Naylor *et al.* 2009), given correlations in the data and Poisson errors in counting integer frequencies associated with a finite temporal sample (Greenhough & Main 2008). In the absence of a good statistical fit to time-limited earthquake datasets, an upper bound can be estimated from the tectonic deformation rates as a constraint (Main & Burton 1984). Tectonic deformation rates are, in fact, remarkably stationary, testament to the relatively constant slow forcing owing to plate tectonics operating on scales of millions of years (DeMets 1995).

The power-law component of the gamma distribution holds over a very wide bandwidth of scales, implying large *θ*. Broad-band scale invariance is one of the most important hallmarks of ‘self-organized criticality’ (SOC), as defined by Bak *et al*. (1987), and first applied to earthquake models by Bak & Tang (1989). However, the converse is not true, i.e. power laws can also arise from systems that are not near critical (Reid & Hughes 2002), and commonly arise in turbulence in the inertial range, in a state very far from SOC. SOC is not a precisely defined concept, but other key hallmarks relevant to the earthquake problem (summarized by Main 1996) include a constant forcing (tectonic strain) rate (consistent with palaeomagnetic and modern satellite data; DeMets 1995); a system that has already self-organized to a state of near-critical stress everywhere, not just at mega-faults (consistent with man-made stress perturbations inducing seismicity remote from plate boundaries; Gupta & Chadha 1995); and small fluctuations about the near-critical state (consistent with the relatively low stress drop observed at all magnitudes; Abercrombie & Leary 1993). The main current debate is whether or not SOC (albeit defined in this loose way) is a stationary state with no predictability of individual events, or there remain systematic fluctuations about this state in an ‘intermittently critical’ way that could lead to the possibility of forecasting individual events in advance, above that expected by say a temporally random process with aftershocks or more generally triggered events. Here, we address a different question—what might determine how close the system average is to the strictly defined critical point, and hence how predictable might we expect the system to be as it fluctuates around its average state?

For at least some non-equilibrium systems, the maximum entropy solution occurs where entropy *production* is also an extremum. For example, in a linear, near-equilibrium system, such as heat conduction through a solid with fixed temperatures at its boundaries and a few degrees of freedom, it can be shown that the entropy-production rate is a minimum (Nicolis & Prigogine 1989; Kleidon & Lorenz 2005). In contrast, when heat is transported by convection in a fluid, the effective heat conductivity of the medium becomes a variable that depends on the fluid velocity, which, in turn, depends on a temperature gradient that is no longer fixed. This inherent nonlinearity in a system maintained in a steady state, far from equilibrium, with feedback, open boundary conditions and many degrees of freedom, can result instead in a steady state of maximum entropy production (MEP) (Kleidon & Lorenz 2005; Whitfield 2005). For example, the temperature profile predicted by MEP is consistent with that of the atmospheres on Earth, Mars and Titan, and perhaps that of mantle convection (Paltridge 1978; Ozawa *et al.* 2003). Critical in this analysis is the ability of the convecting medium to self-tune its effective thermal conductivity.

Dewar (2003, 2005) attempted to generalize the MEP hypothesis to non-thermal systems, deriving a specific generic formulation for information entropy production that includes equivalent flux (analogue thermal dissipation) and local source (frictional dissipation) terms to the thermodynamic equivalents, and is a corollary of the maximum entropy principle under a given set of constraints. Using this theoretical approach, he was the first to connect a maximum in information entropy production to SOC. For example, in the classical sand-pile model (Bak *et al*. 1987), he showed the MEP solution for a constant slow forcing rate (‘grain’ rate) implied a divergence of the scale of fluctuations (equivalent to *θ* defined above), one of the key hallmarks of SOC.

In fact, it has long been suggested that MEP applies to many physical, chemical and biological systems under specific conditions (reviewed by Martyushev & Seleznev 2006). MEP is not a generic probabilistic solution, in the sense that it is always subordinate to stronger macroscopic constraints—it is not a ‘principle’ as such. Explicitly, Martyushev & Seleznev (2006) state that ‘attempts to derive the maximum entropy principle have so far been unconvincing because they often involve additional hypotheses’. Frederiksen & O’Kane (2008) suggest that this is not surprising, as far-from-equilibrium phenomena are often characterized by multiple stationary states and hysteresis. Martyushev & Seleznev (2006) also explain the apparent paradox that minimum entropy production can be viewed as the maximum *possible* entropy production rate for a system otherwise strongly constrained to be near equilibrium, for example, heat conduction in a solid. They provide a historical perspective that provides examples of earlier and alternate mathematical treatments to that of Dewar, and highlight the connection between MEP and SOC as a promising current line of inquiry.

Before exploring this, it is useful to explicitly restate what the concept of SOC as coined by Bak *et al*. (1987) involves. Their idea was that a threshold system, driven far from equilibrium by a constant flux, could spontaneously organize themselves (without tuning) to a ‘critical’ state, i.e. always near the threshold for system-sized failure, characterized by an order parameter (e.g. change in the angle of repose in an avalanche or stress drop in an earthquake) near zero, a diverging correlation length (e.g. in the avalanche or rupture area) and associated broad-band scale invariance leading to power-law scaling statistics (Main 1996). Finite-order parameters or correlation lengths allowed in this slightly loose definition of SOC may be due to finite-sized effects (in which case, the system really is at a thermodynamic critical point) or to near but strictly subcritical behaviour.

Main & Naylor (2008) previously explored the connection between MEP and SOC in model and natural seismicity. They defined entropy production in terms of the ratio of strain energy flux and a ‘temperature term’ not related to heat (particle velocity), but instead to fluctuations in the local strain energy field on a fixed lattice (Rundle *et al.* 1995; Main *et al.* 2000; Sornette 2006, §7.4). Main & Naylor (2008) then applied this to the non-conservative Olami–Feder–Christensen (OFC) (Olami *et al*. 1992) numerical model for earthquakes, which includes an additional degree of freedom in terms of local stress dissipation, and used it to show that MEP occurs strictly in a state of self-organized subcriticality, consistent with all of the key aspects of natural seismicity quoted above, as well as the general observation *b*≈1.

In this paper, we repeat the analysis of Main & Naylor (2008), but using Dewar’s (2003, 2005) broader definition for entropy production applied to a non-conservative system. We first show that, at steady state, the flux and source terms for entropy production in the OFC model are, in fact, equal, and hence that the conclusions of Main & Naylor (2008) are robust with respect to this difference in definition. We then show how the degree of dissipation affects the spatial, temporal and rupture frequency-size distributions in the transition from MEP to minimum entropy production as the conservation parameter is increased to its critical value of unity in the OFC model. The precise critical point, defined for an infinite system by diverging rupture area (or correlation length) and zero seismic efficiency (a dimensionless order parameter defined by the ratio of stress drop to average stress), is in a state of minimum entropy production, associated with a high degree of spatial disorder, a memory-less strain energy field with no predictability and Boltzmann temporal fluctuations in the macroscopic strain energy field. MEP, by contrast, occurs in a near but strictly subcritical state, when there is some memory of past events in the form of correlated ‘domains’, a component of quasi-periodic ‘saw-tooth’ fluctuations in macroscopic strain, low but finite seismic efficiency or stress drop and broad-band scale invariance with an exponent *b* near the observed value of 1, all as seen in natural seismicity.

## 2. The Olami–Feder–Christensen model

The OFC model is a square lattice model of *q* elements, where each cell *i* represents a block connected to its nearest neighbours by a connecting spring with a spring constant *K*_{C} (Olami *et al*. 1992, fig. 1). These blocks rest on a flat plate and are driven at constant strain rate through leaf springs of length *l*_{ 0} and stiffness *K*_{L} that connect them to an upper plate. The two plates represent the sides of the fault zone. Each cell is initialized with a random scalar strain *ε*_{i}. Strain is then accumulated at a constant velocity from the driving plate until a cell reaches its failure stress , at which point the cell ruptures, resets its strain and redistributes a proportion of its strain Δ*ε*=*α**ε*_{i}=*β**ε*_{i}/4 to its four nearest neighbours, where *α*=*K*_{C}/(*K*_{L}+4*K*_{C}) is the elastic parameter of the OFC model and *β*=4*α* is a conservation factor (of local strain). For complete conservation at finite *K*_{C},*β*=1;*K*_{L}=0. If any of these neighbouring cells are now above the threshold, they too rupture in the same manner until no cells are above the threshold. We allow only one event to nucleate in a given time step, equivalent to a condition of slow forcing rate. In accordance with the original OFC model, cells that have been ruptured can accumulate stress from subsequent updates in the same rupture propagation event, i.e. a cell, once ruptured and reset, heals immediately. This introduces a stochastic roughening of the strain field behind the rupture front that is a source of strong spatial fluctuations that increasingly ‘wipe’ the memory of the rupture as a coherent domain as *β* increases (Naylor & Main 2008). We use a spatially uniform fault with for all *q*=200×200 cells, and assign a local reset value of *ε*_{i}=0.

## 3. Entropy production rate

Using the constraints of unit probability, and assuming a stationary state (averaged over the long term) of constant internal energy and mass, and constant energy and mass flux applied to the boundary of a system, Dewar (2003) used the maximum entropy technique applied to path integrals for a Hamiltonian system to derive a general formula for information entropy production in an open, dissipative, non-equilibrium, steady-state system. For a system with a constant boundary flux *F*, a forcing gradient ∇*ϑ* and an internal source term *Q* at inverse temperature *ϑ*, the entropy-production rate is
3.1
The first term describes the rate of the entropy production for the system’s fluxes, and can be interpreted as a product of input fluxes and their conjugate forces. The second term relates to the entropy production owing to frictional dissipation. For the sand pile, the relevant input flux term *F* is the mass flux of sand particles added to the pile (grain rate), and the forcing gradient ∇*ϑ* is proportional to the slope of the sand pile (Dewar 2005).

Assuming a single, stable long-term average state, and using the same constraints as those to define equation (3.1), Dewar (2003) then shows that the maximum probability path, taking into account both reversible and irreversible components, is one of MEP (implicitly subject, as Martyushev & Seleznev (2006) noted later, to additional constraints when they are present). Dewar (2003, 2005) also showed that, in the limit of low driving rate, MEP is associated with an average sand-pile slope near the critical angle of repose and divergence in the variance of the magnitude of the output flux. Both of these are characteristic signals of SOC, the latter with occasional system-size avalanche events. The derivation of MEP is such that it applies only to the statistical properties of the (long-term averaged) stationary state, and not to the individual fluctuations we may otherwise wish to predict. It is therefore more useful in our example in addressing questions such as ‘what is the probabilistic seismic hazard based on stationary plate tectonics’ rather than ‘when is the next big earthquake’.

In the two-dimensional version of the OFC model, the input flux is determined by the constant strain rate *d**ε*/*d**t* from the driving velocity applied to an upper plate. Thus, we can replace the forcing term *F* by a term proportional to *d**ε*/*d**t* (analogous to the mass flux in the sand pile). The most appropriate conjugate gradient term ∇*ϑ* is then proportional to the average strain on the leaf springs 〈*ε*〉 (analogous to the slope in the sand pile). The integration over the volume *V* in equation (3.1) is replaced by a sum over the area *A* represented by the *q* elementary blocks on the fault plane in the OFC model. The forcing term (strain energy flux averaged over many cycles) is then .

The energy associated with local sources of frictional dissipation (earthquakes) is
3.2
where Δ*N*/Δ*t* is the seismic event rate, 〈*E*_{S}〉 is the mean radiated energy per event for the earthquake population and *η* is the seismic efficiency, defined by *E*_{S}=*η**E*, where *E* is the average strain energy held in the crust over time. For elastic radiation, the seismic efficiency is related to stress drop Δ*σ* and the mean stress (average of stresses at and after failure) 〈*σ*〉 by *η*=0.5Δ*σ*/〈*σ*〉. Other constraints may be useful in geodynamic problems in general, for example, in heat flux, momentum, angular momentum, enstrophy etc., but these are not present in the linear elastic, massless OFC model and are not possible to analyse using the recorded properties of earthquakes addressed here, which we show can be explained solely from considering the energy flux and internal source terms.

We quantify an appropriate inverse temperature term *ϑ* in terms of local fluctuations in the strain energy field, determined, in turn, by local differences in the leaf spring strain *δ**ε*. There are twice as many connecting springs as block elements in the two-dimensional OFC model, so the temperature term representing the local fluctuations in strain energy on the fault plane, averaged over the *q* cells in the cellular automaton, is then
3.3
The driving plate is assumed rigid, so that its internal stiffness is infinite, , and hence *ϑ*_{plate}=0. The gradient term is then ∇*ϑ*=(*ϑ*−*ϑ*_{plate})/*l*_{0}=*ϑ*/*l*_{0}.

Including a Boltzmann-like term *k* to retain appropriate dimensions, entropy production is then
3.4
At steady state, Main & Naylor (2008) showed that
3.5
Empirically, the driving stress is related to strain rate in Earth materials under semibrittle conditions across a wide range of scales by a power law *d**ε*/*d**t*=*C**σ*^{n}=*D*(*K*_{L}〈*ε*〉)^{n}, where *C*, *D* and the exponent *n* are constants for a given medium, with *n* estimated to be in the range 2–6 (Carter & Kirby 1978; Newman & White 1997). Applying this to equations (3.4) and (3.5) yields
3.6
Apart from the factor 2, this is identical to equation (3.5) in Main & Naylor (2008), reflecting the fact that, for this problem at steady state, . It follows that all of the main results of Main & Naylor (2008) are robust with respect to this change. We now explore these in more detail.

## 4. Maximum entropy production and self-organized (sub)criticality

Figure 1 shows the mean rupture area, seismic efficiency and entropy production as a function of the conservation factor *β*. Apart from finite-size effects (finite *q*), the mean rupture area diverges and the seismic efficiency tends to zero as *β*→1. The mean rupture area is potted as an analogue for the correlation length of rupture. The seismic efficiency is by definition a dimensionless measure of the difference in energy between ‘broken’ and ‘intact’ phases. It is therefore an appropriate ‘order parameter’ (Stanley 1971), analogous to the density contrast between liquid and vapour phases at the critical point. Diverging correlation length and a zero-order parameter (apart from finite-size effects associated with the model boundary) define a precise critical point at *β*=1 in figure 1.

We find entropy production is a maximum for 0.60<*β*<0.85, depending on the value of the nonlinear exponent *n* (figure 1*c*), with more ‘brittle’ behaviour (higher *n*) corresponding to a more critical system. Thus, MEP occurs in a near but strictly subcritical state, with large but finite correlation length and a low but finite-order parameter near the critical point. Entropy production is a local minimum at *β*=0, and an absolute minimum at the critical point *β*=1. From equation (3.6), the zero entropy production at *β*=1 is a consequence of the extreme case *K*_{L}=0 (equivalent to zero rigidity modulus and *η*→0 in an infinite system). This property of the conservation factor (for which there is no equivalent in the conservative sand-pile model) explains why MEP occurs slightly below the critical point for the OFC earthquake model, rather than at the critical point in the sand-pile model (Dewar 2005). Here, the attractor is strictly in a state of self-organised (sub)criticality.

## 5. Memory and a potential mechanism for characteristic earthquakes

Figure 2 shows how increasing the conservation factor results in a steady decline in the memory of past events in the transition from MEP to minimum entropy production as *β*→1. The memory is expressed as long-lived correlated domains of cells with equal strain, separated by ‘domain walls’ with much stronger local strain gradient where most of the strain energy in the connecting springs is contained at relatively low *ϑ*^{−1}. This solution is energetically favourable because the self-organization at low *ϑ*^{−1} minimizes the macroscopic energy held in the connecting springs by self-tuning to have a large proportion of neighbouring blocks with the same strain. This has direct analogues in the Ising model for magnetism below the Curie temperature (Bruce & Wallace 1989, fig. 8.2).

For intermediate *β*, there is a strong memory of previous large ruptures, resulting in a tendency for a similar large rupture to occur in the future as the load is increased. This is consistent with the hypothesis of spatially ‘characteristic’ earthquakes, where similar extreme events are thought to repeat because of persistent asperities or barriers along the fault plane, analogous to the domain walls of high strain gradient in figure 2. In this sense, the future event ‘knows’ how big it is going to be once nucleated. As the conservation factor and temperature term increase, the spatial field becomes increasingly random owing to the re-rupturing mechanism, until the spatial memory completely disappears at the critical point *β*=1.

Figure 3 shows the time series for the macroscopic strain for the same values of *β* as in figure 2. The spatial memory also results in quasi-periodic saw-tooth fluctuations in the macroscopic strain that gradually disappear as *β*→1, a property first noted by Janosi & Keresz (1993). The combination of repetition of rupture on a persistent large fault segment and quasi-periodic recurrence was used to define the concept of a characteristic earthquake (Schwartz & Coppersmith 1984), based on comparison of earthquake recurrence estimated from catalogues and on longer time scales by geological estimates of repeated slip of similar amounts on exposed fault traces. At the critical point, the fluctuations in macroscopic strain that determine the seismic efficiency also diminish in amplitude and become more Boltzmann-like (figure 3*b*), indicating that the temporal, as well as the spatial, predictability diminish as *β*→1.

In analysis of real data, the concept of characteristic earthquakes remains controversial however. The recent Parkfield earthquake in California was an acknowledged failure of the hypothesis in prospective mode (Bakun *et al.* 2005), and a larger test based on forecasts made in 1980 failed to reject the null hypothesis of spatially localized, but temporally random, occurrence (Kagan & Jackson 1991). Another aspect of the characteristic earthquake model is the elevated occurrence of large events compared with the Gutenberg–Richter trend at low and intermediate magnitudes (Schwartz & Coppersmith 1984). Recently, Naylor *et al.* (2009) demonstrated that very large samples would be required to reject the null hypothesis of a Gutenberg–Richter frequency-magnitude law, in favour of the characteristic earthquake model, owing to integer counting errors for finite temporal samples. Thus, while MEP provides a possible theoretical mechanism for characteristic earthquakes and a degree of predictability, in practice, this may be hard to establish in natural data owing to the proximity to criticality, a lack of direct observation of the stress field with the kind of resolution shown in figure 2 and to sampling uncertainties. For further discussion of the implications of SOC for earthquake predictability, see the Nature website debate at www.nature.com/nature/debates/earthquake/equake_frameset.html.

## 6. Distribution of fluctuations: rupture area

Figure 4 shows the frequency–rupture area statistics for different values of the conservation factor. When 0.60<*β*<0.85, it has previously been reported that the Gutenberg–Richter *b* value is relatively stable near the observed value *b*≈1 (Lise & Paczuski 2001) when the model is driven to steady state. This is confirmed by figure 4, with *b* closest to unity when *β*=0.7 and 0.8. The *b* value drops sharply thereafter, as the largest events reduce the space available for smaller events. A strong finite-size effect is introduced as *β*→1, manifested by curvature in the frequency statistics for the large rupture area (Boulter & Miller 2003). Note the increasing scatter owing to the increasing ‘counting error’ for the rarer largest events, consistent with the analysis by Greenhough & Main (2008).

## 7. Discussion

The ability of a system to achieve any extremum principle is dependent on there being sufficient degrees of freedom for the system to exploit in exploring the phase space. In this paper, the analogue to the thermal conductivity term as a degree of freedom is the conservation factor *β*, which, in turn, depends on the spring stiffness ratio. For . In this model, the connecting springs have the same stiffness parallel and perpendicular to the slip direction. In the slip direction, the connecting springs deform only axially in the conceptual model and the leaf springs only in shear. Hence, the ratio of spring stiffness in the driving direction (connecting to the leaf spring ratio) is equivalent to the ratio of axial modulus on the fault to shear modulus in the surrounding medium. For a Poisson solid, the ratio of axial to shear modulus is around three, equivalent to a compressional to shear wave velocity ratio of √3, a Poisson’s ratio of 0.25 and *β*=12/13, i.e. a value closer to the critical point than the MEP solution in figure 1 for the range of *n* shown. For a given rock type, seismic velocity and Poisson’s ratio vary significantly with crack density (Walsh 1965; Kachanov 1986). The crack density is much higher on or near the fault than around it (Shipton *et al.* 2006), leading to velocity contrasts of 5–15% at typical earthquake nucleation depths (Ben-Zion *et al.* 2007). This is consistent with a softening of *K*_{C} relative to *K*_{L} (of the order of 10–32%), providing a mechanism for reducing *β* in the approach to the MEP solution, and moving the system away from the strict critical point.

A big limitation of the OFC model is to assume the fault plane has already localized to a two-dimensional plane. In future work, it would be interesting to see how softening associated with damage may also be a mechanism for MEP as a driver for deformation localizing from an initially intact body into a macroscopic fault plane (e.g. as in the modified fuse network model of Cowie *et al.* 1993).

## 8. Conclusion

We have derived an expression for the entropy-production rate for the OFC model for earthquake populations using Dewar’s formulation, including both flux and source terms. At steady state, the two are equal for the OFC model, so previous conclusions using only the source term remain robust. Assuming the commonly observed power-law feedback between remote boundary stress and strain rate at steady state, the model maximizes entropy production in a near but strictly subcritical state, with a low but finite seismic efficiency (order parameter), an upper magnitude cut-off (related to the correlation length of ruptured domains) that is large but finite. The diminishing order parameter and diverging correlation length formally identify a critical point where the conservation factor *β*=1. The MEP solution exhibits the universally observed Gutenberg–Richter *b* value of around 1 in frequency-magnitude data. In this state, the model stress field organizes into coherent domains, providing a physical mechanism for retaining a finite memory of past events, with quasi-periodic saw-tooth fluctuations in the macroscopic strain, implying a finite degree of predictability, albeit strongly limited theoretically by the proximity to criticality and practically by the difficulty of directly observing the Earth’s stress field at an equivalent resolution.

## Acknowledgements

M.N. was funded by an EPSRC grant GR/T11753/01 as part of the NANIA project, http://www.ph.ed.ac.uk/nania/. I.G.M. and M.N. benefited from numerous discussions on complexity science with the NANIA group, notably Graeme Ackland and Richard Blythe, who also commented on an earlier draft, and Michael Cates. We also thank Roderick Dewar for comments on a previous version of the manuscript, Alison Ord for inviting us to contribute to the special issue and two anonymous reviewers for their concise and insightful critiques.

## Footnotes

One contribution of 17 to a Theme Issue ‘Patterns in our planet: applications of multi-scale non-equilibrium thermodynamics to Earth-system science’.

- © 2010 The Royal Society