A Carbon-Cycle Based Stochastic Cellular Automata Climate Model
AA CARBON-CYCLE BASED STOCHASTICCELLULAR AUTOMATA CLIMATE MODEL
KLAUS LICHTENEGGER,
Institut f¨ur Physik, Karl-Franzens-Universit¨at Graz,Institut f¨ur Analysis und Computational Number Theory,Technische Universit¨at Graz; A-8010 Graz, Austria,[email protected]
WILHELM SCHAPPACHER,
Institute of Mathematics and Scientific Computing, Heinrichstraße 36Graz University, 8010 Graz, [email protected]
March 23 rd , 2011 Abstract
In this article a stochastic cellular automata model is examined, whichhas been developed to study a “small” world, where local changes maynoticeably alter global characteristics. This is applied to a climate model,where global temperature is determined by an interplay between atmo-spheric carbon dioxide and carbon stored by plant life. The latter can berelased by forest fires, giving rise to significant changes of global conditionswithin short time.
Keywords : Cellular Automata; Forest Fire; Climate Dynamics; CarbonCycle
PACS Nos. : 64.60.ah, 89.75.-k, 89.75.Fb, 89.75.kd
Cellular Automata (CA), as studied first by von Neumann and Ulam [1], providesystems both interesting from a fundamental point of view and for practicalpurposes of model builders. Among the systems most frequently studied areCA versions of daisyworld [2, 3, 4] and forest fire models [5, 6, 7, 8, 9].The main idea is to study a grid or network of identical units (cells), whicheach can take discrete states, and where the time evolution of each cell is gov-erned by simple local rules: The state s i ( t ) of a cell labelled by i at timestep t a r X i v : . [ n li n . C G ] M a r s only determined by its previous state s i ( t −
1) and the states s j k ( t −
1) of asmall number k of neighbours, labelled by j k .There are different ways how to precisely define a neighbourhood, but typ-ically the qualitative behaviour is typically unaffected by the details [10, 11].This can be seen as a footprint of universality – the fact that microscopic de-tails do not influence features of systems which are “critical” in the sense ofWilson [12]. The fact that all rules are local implies a finite maximum speed forthe propagation of all effects.Originally the rules governing cellular automata were strictly deterministic.One can extend this formalism (and make it both more flexible and more real-istic) when including stochastic elements, i.e. give probabilities how the stateof a cell depends on its own previous state and those of its neighbours. On canalso replace the few discrete states by a (quasi)continuum to gain additionalflexibility. The system studied in this article belongs to this class of Stochas-tic Cellular Automata (SCA) models, which the authors apply to a particularfeature of climate-relevant dynamics. It is well-known that the global temperature depends strongly on the amountof greenhouse gases contained in the atmosphere. While H O vapor is themost important of them, its dynamic is extremely complicated and hard totreat reliably in models. Since the amount of H O vapor strongly depends ontemperature, the water vapor cycle may largely act as a positive feedback cycle,enhancing the effects of other greenhouse gases.For these reasons the authors focus on CO , being the most important othergreenhouse gas. Atmospheric CO is crucial for plant life and one typicallyexpects enhanced growth for larger amount of CO available. This is indeedsometimes used as an argument against global warming [13].But in particular when accompanied by climate changes, enhanced plantgrowth can also mean increased vulnerability to large-scale forest fire, which canrelease considerable amounts of CO in very short time [14] (though balancingaspects seem possible [15]).There exists a number of realistic climate models, which attempt to directlydescribe the terrestric system [16]. These models, which are by now quite suc-cessful, are mostly based on the (simplified) treatment of atmosphere and oceandynamics, but mostly neglect the plant life Carbon cyle. Therefore the presentmodel, though significantly more simplistic, can be regarded as complementaryto these approaches. 2 The Model
In this model one studies patches of land, which are, in the course of time,inhabited by plant life. During their growth process, these plants can storecarbon from the atmosphere. The growth rate depends both on temperature(determined by geographical latitude and carbon dioxide in the atmosphere) andthe amount of carbon available (i.e. again carbon dioxide in the atmosphere).For each cell, optimum growth temperature and temperature tolerance arefirst set randomly and then subject to a simple evolutionary process. Uninhab-ited cells can become inhabited by spreading from neighbours. Their character-istic parameters are determined as a weighted average of those of all neighbours,affected by mutations.In addition, there are spontaneous outbreaks of forest fires, and patches aremore vulnerable to them, if they are more heavily overgrown. The fires canspread, and burning cells immediately release all carbon stored in them to theatmosphere in form of carbon dioxide.
Remark : Since the amount of oxygen in the atmosphere is not noticablychanged during all processes, in the following the distinction between “carbon”and “carbon dioxide” will be dismissed and always the word “carbon” is used,where it is implicitly assumed that a sufficient amount of oxygen is availablefor all processes and that carbon is deposited in the atmosphere in the form ofcarbon dioxide.In the process of biological degradation and during fires, not only CO , butalso considerable amounts of H O are released. While, for reasons discussed insection 1.2 we do not attempt to explicitly include water vapor effects, they maybe – at least at some timescales – partially covered too.
The system studied in this approch is closely related to the Clar-Drossel-Schwablforest-fire model [7, 8] which also combines slow growth with rapid burning ofclusters. In contrast to the Clar-Drossel-Schwabl model, the present approachfeatures stochastic propagation of the fire and – more important – the couplingto a global variable, to be interpreted as the atmospheric carbon dioxide contentof the model world.
The model is set up on a rectangular grid, described by several [ N i × N j ]-matrices, which are summarized and explained in table 1. Each cell represents apatch of land, possibly inhabited by a certain amount of plant life with individualcharacteristics. Boundary conditions are chosen periodic in the first ( i ) and openin the second ( j ) direction. 3atrix range meaning L (cid:96) ij ∈ { , } life/inhabitance Matrix C c ij ∈ [0 ,
1] amount of carbon stored in a specific cell T τ ij ∈ R optimal temperature for a cell S σ ij ∈ R + temperature tolerance (width of distribution) B b ij ∈ { , } flag for burning cell (used only in forest fire phase)Figure 1: System matricesA simulation run is performed as follows: Initialization : Cells are inhabited with starting density ρ start , the parame-ters of inhabited cells are determined as c ij = U , τ ij = U and σ ij = 1+ µ · (2 U − U taken individually from a uniform [0 , N steps times:1. Spreading : Each cell with at least one inhabited neighbour has a chanceto become inhabited as well. The probability is given by p ( i, j ) = 14 + 4 π M (cid:88) i (cid:48) j (cid:48) (cid:48) c i (cid:48) j (cid:48) . (1)where the prime at the sum indicates summation over next neighbours,with diagonal elements ( i ± , j ±
1) weighted with the Moore factor [10] π M = . In case the cell becomes inhabited, its parameters are determinedas τ ij = 1 N (cid:88) i (cid:48) j (cid:48) (cid:48) c i (cid:48) j (cid:48) τ i (cid:48) j (cid:48) + µ · (2 U − , (2) σ ij = 1 N (cid:88) i (cid:48) j (cid:48) (cid:48) c i (cid:48) j (cid:48) τ i (cid:48) j (cid:48) (1 + µ · (2 U − , (3)with N := (cid:80) (cid:48) i (cid:48) j (cid:48) c i (cid:48) j (cid:48) and U individually taken from a uniform randomdistribution on [0 , τ j = 1 (cid:80) i (cid:96) ij N i (cid:88) i =1 (cid:96) ij τ ij and σ j = 1 (cid:80) i (cid:96) ij N i (cid:88) i =1 (cid:96) ij σ ij (4)are stored for later time series analysis. In principle, one could introduce an additional multiplicative parameter α spread ∈ (0 ,
1] infront of the sum, the same way it is done for the following processes. Our choice of α spread = 1effectively fixes the timescale. Accordingly the rates of all other processes (except spreadingof forest fires) are determined relative to this spreading rate. Determining global variables (1) : The carbon content C of the atmo-sphere is determined as C = 1 − N i N j (cid:88) i,j c ij (5)and the temperature is determined as T j = α Tvar j − N j − − α Tvar ) C . (6)The value of C is stored for later time-series analysis.3. Growth : Each cell growths with its individual rate, c ij ( t + 1) = med (cid:110) , c ij ( t ) + α growth C ϕ ( T j , τ ij , p ij ) , (cid:111) (7)(where “med” denotes the median value) with the temperature adaptationmodel function ϕ ( T j , τ ij , σ ij ) = N exp (cid:26) − ( T j − τ ij ) σ ij (cid:27) − h hostil , (8)where the normalization N = N ( τ ij , σ ij ) is chosen such that (cid:90) ϕ ( T j , τ ij , σ ij ) = 1 − h hostil . (9)This simple Gaussian shape already shows the most important featuresone would expect of such a function. It allows both extreme specializationto a certain temperature and adaptation to a wider range. For h hostil > c ij ( t + 1) to be = 0, the cell dies.4. Determining global variables (2) : Again Carbon content and temper-ature are determined according to equations (5) and (6).5.
Forest fires : Each inhabited cell ( i, j ) has a chance to catch fire, withprobability p = α fire T j c ij . (10)If there any burning cells (cells that caught fire), a complementary simu-lation on a short timescale is started. The fire may spread from burningcells ( b ij = 1) to their neighbours:(a) In each step of the burning process, each inhabited cell can catch firewith probability p = T j c ij
14 + 4 π M (cid:88) i (cid:48) j (cid:48) (cid:48) b i (cid:48) j (cid:48) . (11)5aram. range meaning N i ∈ N ≥ linear size of the grid in first direction (open b.c.) N j ∈ N ≥ linear size of the grid in second direction (periodic b.c.) t max ∈ N simulation time ρ start ∈ (0 ,
1] initial density of inhabited cells α fire ∈ R + probability of a cell spontaneously catching fire α growth ∈ R + measure for growth rate h hostil ∈ R ≥ “hostility” of environment α Tvar ∈ [0 ,
1] Geographic variability of the temperature (6) µ ∈ [0 ,
1] global mutation rateFigure 2: A summary of the system parameters, which are explained in detailin sec. 2.3(b) In each individual cell, the fire lasts for one timestep and is extin-guished afterwards ( b ij → The result of a typical simulation run is displayed in figure 3. After a simulationrun, the following data is available: • The time series of carbon C in the atmosphere, • the time series of average temperature adaptation A ( t ) := 1 − N L ( t ) (cid:88) ij (cid:96) ij ( t ) | τ ij ( t ) − T j ( t ) | , (12) • the time series of average versatility V ( t ) := 1 N L ( t ) (cid:88) ij (cid:96) ij ( t ) σ ij ( t ) , (13)where the defintion N L ( t ) := (cid:80) ij (cid:96) ij ( t ) has been used. From these time series,the authors have extracted maximum, average and minimum values for A and6igure 3: Results of a simulation run for the parameter choice N i = 12, N j = 24, t max = 256, ρ start = 0 . α fire = 0 . α growth = 0 . h hostil = 1, α Tvar = 0 . µ = 0 . τ -values, the current σ -values and a timeseries of average temperature adaptation (12).The fourth row displays the time series of row-averaged τ and σ and a timeseries of average versatility (13). (In case that row-averaged values were notwell-defined for certain times [when no cell of that particular row was inhabited],this is reflected by white spots in the corresponding graphs, cf. figure. 8.)7 and certain characteristics of the Fourier spectrum of C (examining onlyabsolute values, i.e. the amplitude spectrum X ( n ), which is closely related tothe power spectrum).These characteristics are the position and the value of the dominant maxi-mum beyond n = 1 and the parameters of a power-law fit P ( n ) ≈ Cn α + o offset . (14)In addition the authors have analyzed the number of exceptional increasesof C (usually caused by forest fires), which are defined as timesteps τ , where C ( τ ) > C ( τ − C ( τ − < C ( τ − C ( τ ) := C ( τ ) − C ( τ − Comparison of simulation runs for various sets of parameters revealed that themost crucial parameters are α fire and h hostil . Thus the authors have performedsimulations for values of these two parameters chosen on a grid. For a widemesh size and α Tvar = , the results are shown in figure 4.Some of the most important features of the model are already clearly recog-nizable: For α fire = 0 there are no forest fires, so for small h hostil all carbon isstored in the cells. For large values of h hostil , it is still possible that the wholepopulation dies out due to lack of adaptation. This happens in general (for thechosen set of parameters) for h hostil ≥ .
75, regardless of α fire .For h hostil = 0 and α fire small, but nonzero, one finds a local maximumof atmospheric carbon. This peak is a key feature of the model and will beexamined closely in the following analysis. The basic explanation is that verysmall α fire implies long periods without fires, in which large amounts of carbonare stored, such that one forest fire can have devastating consequences andeliminate almost all cells at once, releasing large amounts of carbon.The authors find adaptation maxima at h hostil ≈ h hostil the punishment for bad adaptation is less severe,while for even larger values of h hostil the population more frequently dies outcompletely without having time to adapt.To study these effects more closely, the authors have restricted the grid toa smaller region in parameter space. Performing again several simulation runsgave the results depicted in figure 5. In addition they changed the parameter α Tvar to zero (figure 6) and one (figure 7) and find that the former case (wheretemperature is solely dominated by C ) is significantly different, while the lattercase gives results very similar to those obtained for α Tvar = .8 h hostil α fire X max M pos αA N C inc ∆ C inc h h α f h h α f h h α f h h α f h h α f h h α f C h hostil α fire X max M pos αA N C inc ∆ C inc h h α f h h α f h h α f h h α f h h α f h h α f σ ( C ) h h α f . σ ( P max ) h h α f . σ ( M pos ) h h α f . σ ( α ) h h α f . σ ( A ) h h α f . σ ( N C inc ) h h α f . σ (∆ C inc ) h h α f . Figure 4:
Results of simulation runs (10 runs for each specific combination of pa-rameters) for the parameter choice N i = 12, N j = 24, t max = 256, ρ start = 0 . α growth = 0 . α Tvar = 0 . µ = 0 .
1. The parameters α f = α fire and h h = h hostil havebeen varied as α fire = 0 . k , k = 0 , . . . ,
20 and h hostil = 0 . (cid:96) , (cid:96) = 0 , . . .
8. On theleft hand side, the main plot shows the average carbon content C of the atmosphere.On the right hand side, the first row of graphs gives the amplitude of the dominantFourier mode, X max , its position M pos and the exponent α of the spectral power-lawfit (14). The second row shows the average temperature adaptation ( A ), the totalnumber N C inc of timesteps for which C increased (which corresponds to the numberof significant forest fires) and the average change ∆ C inc of carbon content during suchevents. Plots of the same quantities are given one more time (second graph on theleft hand side, third and forth row on the right hand side), but usually shown froma different viewing position. In the last row the standard deviations σ of the sevenquantities under consideration are given. h hostil α fire X max M pos αA N C inc ∆ C inc h h α f h h α f h h α f h h α f h h α f h h α f C h hostil α fire X max M pos αA N C inc ∆ C inc h h α f h h α f h h α f h h α f h h α f h h α f σ ( C ) h h α f . σ ( P max ) h h α f . σ ( M pos ) h h α f . σ ( α ) h h α f . σ ( A ) h h α f . σ ( N C inc ) h h α f . σ (∆ C inc ) h h α f . Figure 5: Results of a sequence of 10 simulation runs. Parameters are chosenas in figure 4, but only the most interesting region is examined by choosing α fire = 0 . k , k = 1 , . . . ,
20 and h hostil = 0 . (cid:96) , (cid:96) = 0 , . . .
15. Again the graphsare displayed for two different viewing positions.Note that for larger values of α fire the carbon content C is almost constantfor small h hostil and increases monotonically as a function of h hostil , starting at h hostil ≈
1. For small values of α fire and h hostil , there are large fluctuations in C , to be noticed also in the graph for σ ( C ).10 h hostil α fire X max M pos αA N C inc ∆ C inc h h α f h h α f h h α f h h α f h h α f h h α f C h hostil α fire X max M pos αA N C inc ∆ C inc h h α f h h α f h h α f h h α f h h α f h h α f σ ( C ) h h α f . σ ( P max ) h h α f . σ ( M pos ) h h α f . σ ( α ) h h α f . σ ( A ) h h α f . σ ( N C inc ) h h α f . σ (∆ C inc ) h h α f . Figure 6: Results of a sequence of simulation runs. Parameters are chosen asin figure 5 except for α Tvar = 0. Again the main graphs are displayed for twodifferent viewing positions.For this choice of parameters, the temperature does not depend on geographicallatitude, but only on C . Apparently in this situation, for small values of α fire it is more likely that the atmosperical carbon level remains low, i.e. the cellsinhabitance is stable. This is also reflected in in the Fourier amplitude spectrum,where the maximum still has low intensity, indicating the absence of significantoscillations as they can be recognized in figure 3.11 h hostil α fire X max M pos αA N C inc ∆ C inc h h α f h h α f h h α f h h α f h h α f h h α f C h h α fire X max M pos αA N C inc ∆ C inc h h α f h h α f h h α f h h α f h h α f h h α f σ ( C ) h h α f . σ ( P max ) h h α f . σ ( M pos ) h h α f . σ ( α ) h h α f . σ ( A ) h h α f . σ ( N C inc ) h h α f . σ (∆ C inc ) h h α f . Figure 7:
Results of a sequence of simulation runs. Parameters are chosen as infigure 5 except for α Tvar = 1. Again the graphs are displayed for two different viewingpositions. In this case, where temperature is independent of the carbon content C ,one finds results not so different from those displayed in figure 5. This, togetherwith the significantly different results shown in figure 6 seems to indicate that thesystem characteristics change signifcantly for small values of α Tvar , i.e. the mostly C -dominated case. .3 Interplay of parameters The model parameters are not completely indepent. For example it is clear thatin figure 4 the population typically dies out for h hostil ≥ .
75. One could expect,however, that a larger mutation rate µ may be beneficial in dealing with hostileenvironments and that one could find values of µ that allow the population tosurvive in most cases.This is indeed the case, as illustrated in figure 8. For h hostil = 1 .
75 and µ = 0 . µ = 0 .
5, thoughin this case it repeatedly comes close to extinction. This is a hint that also inthis model a mutation rate which is too large can threaten the survival of aspecies.
There are several possible extensions of this model, which might be interestingto study: • One could modify the grid by of introduction of uninhabitable spots, dis-tributed randomly or as clusters of some fractal shape. One may also con-sider the growth of such uninhabitable spots (which can be interpreted asspreading of cities and deserts). • In a more complete model one could explicitely consider the influence ofsoil, which can store carbon from deceased plants. Another importantfactor is humidity, which can reduce the rate of forest fires and also hasinfluence on temperature. This could help to make the model more real-istic, but at the same time harder to analyze. • One could dismiss the normalization condition C ≤ • One could try to explicitely include the effects of water vapor by setting upa a separate cycle for H O, coupled to the plant life grid, but endowed witha separate dynamics on a shorter timescale. It might even be worthwhileto think about coupling a Stochastic Cellular Automata Carbon climatemodel as proposed in this article to a more traditional model based onatmosphere and ocean dynamics.13igure 8: Large values for h hostil can be partially compensated by increasing themutation rate µ . These plots show typical simulation runs for h hostil = 1 .
75 and µ = 0 . µ = 0 . µ = 0 .
5. (Other parameters were chosen as in figure 3.)14
A fascinating extension may be the coupling of the system to an “ocean”,i.e. a carbon sink with limited absorption rate and temperature-dependentcapacity.This may be especially important, since the solubility of CO in waterdecreases with increasing temperature . So while the oceans can storeenormous amounts of CO , their capacity decreases with global warming.This could lead to a dramatic positive feedback chain:In such a scenario, an increase of temperatures above a certain criticallevel could lead to a massive release of CO from the oceans and thus afurther rise of temperature. Confronted with this frightening perspective,even simple models which allow to study the dynamics of such catastrophicbehaviour can be extremely valuable.Apart from these extensions, there are also several aspects of the original modelthat should be studied in greater depth. The interplay of model parameters, assketched in subsection 3.3, certainly deserves closer examination. In addition itmay be interesting to study size-dependence, in particular the correspondencebetween “critical” α fire and the extension of the grid (for example by consideringthe limit N i N j → ∞ , α fire → N i N j α fire held fixed). The authors have studied a “small-world” model designed to describe the in-terplay between global climate, carbon storage and sudden outbreaks of forestfires.They have seen that the least predictable behaviour is found for small, butnonzero values for the probability α fire of a cell catching fire. While in most casesan increased hostility of the environment leads to a smaller average populationand thus to a higher carbon content of the atmosphere, this is not necessarilytrue in the most complex (and probably most realistic) case of small, but nonzero α fire .This “critical” behaviour is more dramatic in the case of temperature dom-inated by atmospheric carbon (instead of geographic parameters), in terms ofthis model for α Tvar ≈ This is in general true for the solubility of gases in liquids, while it is exactly the otherway round for the solubility of solids in liquids. This is easy to understand, since enhancedthermal movement of the molecules makes it easier for gases to evaporate, while it makes itharder for solids to crystallize. eferences [1] Von Neumann J. et Burks A. ed., Theory of Self-Reproduction Automata ,University of Illinois Press, 1966, p. 77, in Ostolaza J.L., Bergareche A.M.,La vie artificielle, Seuil, Paris, 1997, pp. 37-38. Translated from French.[2] Watson, A.J., J.E. Lovelock (1983).
Biological homeostasis of the globalenvironment: the parable of Daisyworld . Tellus B 35 (4): 286-9. Int. Mete-orological Institute;[3] Booth, Ginger; Classical and CA-based daisyworld simulator, http://gingerbooth.com/courseware/pages/demosdaisy.html [4] Graeme J. Ackland, Michael A. Clark, Timothy M. Lenton,
Catastrophicdesert formation in Daisyworld , Journal of Theoretical Biology 223 (2003)39-44[5] Bak, P., Chen, K. and Tang, C. (1990),
A forest-fire model and somethoughts on turbulence , Phys. Lett. A 147, 297-300.[6] Chen, K., Bak, P. and Jensen, M. H. (1990),
A deterministic critical forest-fire model , Phys. Lett. A 149, 207-210.[7] S. Clar, B. Drossel, F. Schwabl,
Scaling laws and simulation results forthe self-organized critical forest-fire model , arXiv.org:cond-mat/9405008 ,1994[8] S. Clar, B. Drossel, F. Schwabl, Forest fires and other ex-amples of self-organized criticality , COND.MAT. , 6803, arXiv.org:cond-mat/9610201 , 1996[9] R. Krenn, S. Hergarten, Cellular automaton modelling of lightning-inducedand man made forest fires , Natural Hazards and Earth System Sciences, 9(2009) 1743-1748[10] Klaus Lichtenegger,
Stochastic Cellular Automata Modelsin Disease Spreading and Ecology , diploma thesis, 2005; http://physik.uni-graz.at/ ∼ kll/cthesis.pdf [11] Klaus Lichtenegger, Wilhelm Schappacher, Phase Transition in a Stochas-tic Forest Fire Model and Effects of the Definition of Neighbourhood ,IJMPC 20, 8(2009) pp. 1247-1269, arXiv:0902.3680v1 [nlin.CG] [12] K. G. Wilson and J. B. Kogut,
The Renormalization group and the epsilonexpansion , Phys. Rept. (1974) 75.[13] Fowler, Simon, Positive Effects of Carbon Dioxide for Plant Growth , http://ezinearticles.com/?Positive-Effects-of-Carbon-Dioxide-for-Plant-Growth&id=1607 Wildfire Drives Carbon Levels In Northern Forests ; see [15] University of California, Irvine (2006, November 17),
For-est Fires May Lead To Cooling Of Northern Climate , see [16] Randall, D.A., R.A. Wood, S. Bony, R. Colman, T. Fichefet, J. Fyfe, V.Kattsov, A. Piman, J. Shukla, J. Srinivasan, R.J. Stouffer, A. Sumi andK.E. Taylor, 2007: Climate Models and Their Evaluation. In: