Mathematical Modelling of Climate Feedback Printer-friendly Version Interactive Discussion Earth System Dynamics Discussions Mathematical Modelling of Positive Carbon-climate Feedback: Permafrost Lake Methane Emission Case Mathematical Modelling of Climate Feedback Printer-friendly Version Interacti

Introduction Conclusions References Tables Figures Back Close Full Screen / Esc This discussion paper is/has been under review for the journal Earth System Dynamics (ESD). Please refer to the corresponding final paper in ESD if available. Abstract Introduction Conclusions References Tables Figures Back Close Full Screen / Esc Abstract The permafrost methane emission problem is in the focus of attention of different climate models. We present new approach to the permafrost methane emission model-ing. The tundra permafrost lakes is potential source of methane emission. Typically, tundra landscape contains a number of small lakes and warming leads to lake exten-5 sion. We are making use of this process by the nonlinear theory of phase transitions. We find that climate catastrophe possibility depends on a feedback coefficient connecting the methane concentration in atmosphere and temperature, and on the tundra permafrost methane pool.


Introduction
In this paper we consider the methane emission problem, which is important for climate change.Recently it has become that in Siberia permafrost exists a "methane bomb", namely, about 80 × 10 15 kg of carbon in a frozen state.As a result of tundra permafrost thawing, a number of small lakes form and extend, with methane entering the atmosphere.In turn, atmosphere methane can reinforce warming, and there appears a positive feedback that finally can lead to a climate catastrophe.
Firstly, we consider a mathematical model of lake growth and then using this model we propose a simple phenomenological equation that allows us to evaluate the impact of the Siberia permafrost on climate.
Mathematically, permafrost thawing can be described, for example, by the classical Stefan approach (Caginalp, 1989).However, it is difficult to resolve the multidimensional Stefan problem numerically.It is known that two-dimensional Stefan layers are unstable which can lead to difficulties in numerical simulations (Caginalp, 1989).To overcome this difficulty, we can use a modified approach based on the phase transition theory (Caginalp, 1989).This takes into account that thawing layer has a small but non-zero width.Figures

Back Close
Full The main idea is as follows.Typically, the lakes have horizontal sizes of 100-2000 m, and a small depth.It is clear that temperature at the lake bottom is essentially below than that at the lake boundary, where depths are smaller.The 3-D phase transitions are difficult to describe, in general.However, we use some elementary arguments that allows us to obtain simple approximations.
The transition from the frozen state to the thawing state is a microscopic process, while lakes are great macroscopic objects.Thus we can assume that locally a lake boundary is a sphere of a large radius of curvature R. Therefore, the lakes are similar to a shallow large bowl (cup) (see Fig. 1).Moreover, the growth is a slow process.Under such assumptions, thawing front velocity can be investigated.Indeed, there are possible asymptotical approaches based on so-called mean curvature motion (Angenent, 1994;Angenent et al., 1995;De Mottoni and Schatzman, 1995).This approximation can be obtained for all models where we have an energy functional and where the thawing front is narrow.Notice that in this approach the thawing layer has a size, it is not infinitely narrow.The mean curvature motion has been extensively investigated by a great number of papers (Brakke, 1978;Huisken, 1984;Crandall and Lions, 1983;Ilmanen, 1990;Bronsard and Kohn, 1991;Chen, 1992;Barles et al., 1993;Ilmanen, 1993;Buckdahn et al., 2001).In this way we have obtained that the lake increasing rate is proportional to a constant, that can be empirically estimated, plus a small contribution inversely proportional to averaged lake curvature R, the lake surface area A(R) is, approximatively, cR 2 , where c is a constant connected with averaged lake form.According to these theoretical ideas, the lake surface form becomes close to ellipsoidal or circular for large times.
As a result, we obtain a deterministic equation serves as an extremely simplified model of the lake growth.experimental data (Downing et al., 2006).Parameters of this Pareto density can be evaluated by experimental data obtained for the Central Yakutia region (Kravtsova and Bystrova, 2009).At the second step, we use this model to compute the methane emission generated by lakes of the Central Yakutia.As a result, we can, therefore, propose here a simple method to compute methane emission.using a natural assumption that the horizontal dimensions of the lakes are much larger than the lake depth.This allows us to obtain simple relations for the methane increasing rate.
This estimate for the lake increasing allows us to obtain a formula describing a dependence of methane production in tundra on temperature.Here we also take into account empirical relations for methane emission by a lake obtained in Frolking and Grill (1994) and Kotsyurbenko et al. (2004) and some theoretical ideas leading to a relation of Arrhenius's type (Komarov, 2003).
The numerical simulations based on these theoretical results show that there are three different possible scenarios of methane emission.The first case appears if the feedback coefficient γ that connects temperature at the lake surfaces and methane concentration is small enough.Then the methane concentration growth can be described by a linear function in time.
The most interesting situation occurs if γ is large enough.Then one has firstly a linear growth on an initial time interval, then a parabolic increasing curve and finally, this transforms into a sharply increasing curve, that can be interpreted as a "Arctic Armageddon" (Kerr, 2010).For real parameters, however, Armageddon proceeds within a very large time period (years).Introduction

Conclusions References
Tables Figures

Back Close
Full 2 Lake growth model

Stefan problem
Let us consider a cylindrical domain where Ω is a sub-domain of 2-D sphere corresponding to the permafrost surface, (x, y) are coordinates on Ω, z ∈ [0, h].
The unknown function is temperature u = u(x, y, z, t).We assume that h L, where L = diam (Ω) and that the boundary of Ω is a smooth curve.Let us denote by l (x, y, t) the thawing front position, the front surface is defined by z = l .The classical Stefan problem can be then formulated as follows: and we have the following boundary conditions on the depth ũz (x, y, h, t) on surface u(x, y, 0, t) = U(x, y, t), (4) and at the thawing transition front: where n is the unit normal vector with respect to the thawing front, v is the normal front velocity, Q is a dimensionless latent heat.We define a separating surface Γ by the following equation u(x, y, l (x, y, t), t) = θ. (6) Therefore, the water corresponds to the domain {u(x, y, z, t) > θ}, while the ice domain is {u(x, y, z, t) < θ}.Introduction

Conclusions References
Tables Figures

Back Close
Full

Phase transition theory approach
Following Caginalp (1989), one can use an improved approach, which reduces to the Stefan method as the interface width is small.We introduce the so-called order parameter which equals 1 f for the interface.We denote this order parameter by φ, φ(x, y, z, t).We use the following equations for two unknown functions φ, u: where g = 0.5(φ − φ 3 ), θ is the temperature of the phase transition (temperature of thawing).
This model allows us to describe the thawing zone of non-zero width.To explain the main idea of the front formation in this model, let us consider the 1-D variant of Eq. ( 8) removing the term with the temperature: where we assume that the front moves along the z-axis.Let us assume α = 1.Then this equation has the solution tanh(x/ ), where = a −1/2 ξ −1 is a parameter that defines the front width.To describe the influence of 2(u − θ), we can apply the perturbation theory.
The front appears at a point, where u ≈ θ, i.e. at the thawing temperature.
We have the following boundary conditions on the depth ũz (x, y, h, t)

Conclusions References
Tables Figures

Back Close
Full The thawing transition front should be defined now by relation just the initial value problem is resolved.Here α, a, ξ, b are parameters.
Let us consider a connection between this formulation and the Stefan problem (Caginalp, 1989).
Assume the parameters a, ξ → 0, and ξa −1/2 , α are fixed.In this case the Ginzburg-Landau (phase field model) reduces to the modified Stefan problem with where κ is the front curvature, ∆S is a constant describing a difference between termodynamics of ice and water.Additional contributions in the modified and alternative modified situations have a stabilizing influence.
Notice that numerically the phase field approach is more attractive than the method based on the Stefan problem.Indeed, it is difficult to prove existence of solution to the 3-D Stefan problem.The phase field model is well posed and the corresponding initial value problem is well posed.Moreover, here we can take into account that the thawing front has a non-zero width.We will show below that complicated equations lead to quite transparent results, if we take into account actual lake parameters.

A phenomenological model for lake growth
First we make some preliminary mathematical comments.We consider the case C. If ∆s in Eq. ( 15) equals zero, we obtain the relation v(x, y, z, t) = − µκ(x, y, z, t), (16) Introduction

Conclusions References
Tables Figures

Back Close
Full where v is the normal thawing front velocity at the point (x, y, z), κ is the mean front curvature at this point, µ is a positive coefficient.This relation describes so-called mean curvature flows well studied by numerous works (Huisken, 1984;Gage and Hamilton, 1986;Ecker and Huisken, 1989;Bronsard and Kohn, 1991;Ecker and Huisken, 1991;Evans et al., 1992;De Mottoni and Schatzman, 1995;Giga, 2002).Such a flow also occurs if we eliminate the temperature from our equations.Actually, Eq. ( 16) also is a nontrivial nonlinear equation, even in 2-D case.In fact, if the front Γ is defined by equation z = S(x, t) then we have the following formula for κ: Relation Eq. ( 16) entails remarkable consequences that can be checked experimentally.Let us consider the plane case.Then our fronts are curves.All front closed curves, which initially are not too distinguished of ellipses or circles, approaches for large times to circles.For circular fronts of the radius R(t) Eq. ( 16) takes the form and can be resolved explicitly.This solution describes a circle that slowly shrinks.
To take into account a natural lake extension, we add to Eq. ( 18) a positive term and obtain the following asymptotical relation for averaged lake growth: the quantity δ can be expressed via microscopic parameters a, ξ, b, however, it is simpler to find this quantity by experimental data since the δ determines the main contribution in the lake area increasing.This equation can be resolved analytically and solution describes a growing lake of a circular form.Formally, in the case C this constant δ Introduction

Conclusions References
Tables Figures

Back Close
Full can be obtained from Eq. ( 15) if ∆S = 0. To explain δ appearance, we use Eqs.( 20) and ( 21).For small curvatures (large R) we can use an one-dimensional approximation of these equations where u, φ depend on n, t and n(x, y, z, t) is the distance between the observation point M = (x, y, z) and the front.We assume here that n has a sign, i.e. n > 0 if M before the front and n < 0 if M is beyond the front.Then, up to small corrections proportional to O(κ), one has ∇ = ∂/∂n, and the Laplacian ∆ = ∂ 2 /∂n 2 .
Assume, moreover, following standard ideas (Molotkov and Vakulenko, 1988), that the melting front can be described by a slowly moving traveling wave, φ = Φ(Z), where Z = n − V t is a new variable and where the velocity V = δ is small.Then Eq. ( 20), Eq. ( 21) can be transformed into Following some standard calculations (Molotkov and Vakulenko, 1988) one obtains that (u(z)−θ)Φ Z dZ that allows us to find V = δ as a function of the microscopic parameters α, a, ξ, b.
Notice, however, that some of these parameters are difficult to evaluate precisely, it is simpler to find δ directly from experimental data.This constant defines a main contribution to growth of the lake area.

Stochastic effects, Fokker-Planck equation and Pareto distribution
Actually, the lake growth is more complicated process.To take into account possible stochastic effects one can use the Fokker-Planck equation having the form

Conclusions References
Tables Figures

Back Close
Full where ρ = ρ(R, t), R > 0 is a time depending lake radius.To see the form of f , let us remind that for d = 0 the Fokker-Planck equation is equivalent to deterministic equation (Isihara, 1994) dR/dt = f (R).Therefore, in our case f = δ − κ/R.Let us set δ = 0 (the influence of δ > 0 reduces to a trivial shift).It is well known that, typically, ρ(R, t) tends to a stationary equilibrium distribution ρ eq (R) for large times t.
This equilibrium distribution can be defined by the relation This result is consistent with some experimental observations (Downing et al., 2006).
In fact, the lake areas in different regions are distributed according to the Pareto law, where A min is a minimal lake size, k is a positive coefficient.Clearly, A = πR 2 and β = −(2 k + 1).Some experimental data show that the parameter k lies within the interval (0.6, 1).our ideas on the mean curvature motion are correct, then the lakes should have, approximatively, a circular form.This is confirmed by some observations, see for example, aerospace photos from Kirpotin et al. (2008).

Methane emission
To describe the methane emission rate, let us assume, following Komarov (2003), that the probability to emit a gaseous methane molecule is proportional to where u is the temperature, k is the Boltzmann constant, U 0 is a unknown parameter to adjust, which describes a potential barrier for a reaction solid-gas producing methane emission.
Then the total rate of the methane emission per unit time is given by where u = u(x, y, z, t) is the absolute temperature at the point (x, y, z) at the moment t, ρ 0 is the methane density in the domain D M , this domain is the methane emission zone depending on time t, and formally defined by Now the known results of the phase field theory allow us to propose a simple semiheuristic model describing the methane emission by a system of a number of almost plane lakes.We assume that the lakes are shallow and the lake growth reduces to the radius R i increasing.Indeed, the most of such tundra permafrost lakes are shallow

ESDD Introduction Conclusions References
Tables Figures

Back Close
Full (depth 1-10 m, typically 2-3 m) but large (100-1000 m in diameter).The thawing and erosion lead to the diameter extension, the depth increases essentially slower.To describe an averaged lake growth we use Eq. ( 19) assuming that, in average, the form of the lake is close to circular.Indeed, these lakes developed as a result of longtime evolution, thus we can assume (see section on mean curvature) that original, possibly complicated, forms approached to circular.Then, we can use mean curvature motion Eq. ( 19), where µ, δ are positive coefficients.This coefficient can be calculated from microscopic parameters a, ξ, b of the phase transition theory, but it is more reasonable to consider Eq. ( 27) as a phenomenological model, where they should be adjusted by experimental data.
We assume that there are a number of R i , i = 1, 2, ..., N lakes.The averaged lake Let us turn to Eq. ( 19).Then the total methane emission rate per unit time V met is given by This implies, after a simple calculation, that where µ, δ are constants that can be found from observation data, R is the averaged lake radius.According to the previous section, the coefficient C depends on temperature u as Here c 0 , a are new coefficients to adjust (from experiments).Naturally, we should average u over some time period, since the lake processes are slow with respect to atmosphere oscillations, and over a region under consideration.

Conclusions References
Tables Figures

Back Close
Full Finally, we conclude with the following simplest semi-phenomenological relation where β is a cumulative constant that should be found from experimental data, b, d are positive constants, the u is the absolute temperature.This equation is in an accordance with some phenomenological models.We can mention here two following models.The first is a result of observations under a large mire "Bachkarevskoe" within 1994-2002 (Kotsyurbenko et al., 2004).The following empirical formula was proposed (an analogous formula is obtained theoretically by Frolking and Grill, 1994, and a similar relation see Krylova et al., 2001) where F is the flux measured in Mg m −2 per hour, T 12 the temperature at 12 cm depth, W is the water level.To connect Eqs. ( 28) and ( 29) we can use the Taylor extension for the right hand side of Eq. ( 28) at some averaged temperature u = u 0 assuming that the absolute temperature deviations are small with respect to u 0 .This gives us Eq. ( 29).If we take u 0 = 273 K, then b = 0.126 × 273 −2 .Notice that data (Krylova et al., 2001) give coefficient 0.18 instead of 0.126 in Eq. ( 29), i.e. effect is even stronger.Notice that approximation Eq. ( 29) is correct for bounded times, for large times where the feedback is essential, using this relation gives a numerical instability and it is better to use the Arrhenius formula Eq. ( 28).

Numerical results
Here we focuses our attention on a contribution of the Central Yakutia lakes into methane emission and methane concentration growth.Notice that some experimental observations show that many lakes in West Siberia vanish or diminish its surfaces (Kirpotin et al., 2008).The most new data from Kravtsova and Bystrova (2009) and Introduction

Conclusions References
Tables Figures

Back Close
Full  (2011) show that in almost all regions of West Siberia lakes evolve in an extremely sophisticated manner, however, the lake general area conserves, but in the Central Yakutia a strong area growth is observed.A similar picture is observed in some other regions.Authors Kravtsova and Bystrova (2009) explain this phenomenon as a result of higher averaged temperature in these regions.
We simulate a development of a lake system consisting of N 1 lakes.The start area A i (0) were distributed according to the Pareto law, which is adjusted to be consistent with observed parameters.The values A i (t) can be then computed by Eq. ( 19), where we take µ = 0 to simplify computations, and, moreover, for large lakes, this contribution should be small.Then simulations show that within time intervals of order of hundred years the lake area S(t) increases linearly in time, S(t) = S(0) + αt.The parameter α is chosen according to data (Kravtsova and Bystrova, 2009): within 1976-2000 the area is increased from 0.5 × 10 10 m 2 up 1 × 10 10 m 2 .
Moreover, in our simulations we use the following experimental facts (Matthews and Fung, 1987;Cao et al., 1996;Bartlett and Harriss, 1993): the total mass of the methane in atmosphere is about 5 × 10 12 mg.The concentration growth is 1 % per year (Kiehl and Trenberth, 1997).The emission rate is computed by ( 28), and we take into account a small positive feedback assuming that u(t) = ū + γX , where X (t) is the total methane mass in atmosphere, γ is the feedback coefficient, u(t) is temperature, ū is a mean temperature within July, June and August taken without warming effect.
To obtain a rough estimate of the feedback coefficient γ we can use the following arguments.The complete effect of atmospheric greenhouse gases is about 5 K, the methane contribution is within 4-9 %.The second equation takes into account the flow from increasing lake area in the Central Yakutia: where The coefficient 0.25 is inserted since emission goes only within June, July and August.
The simulations give the following typical plots of X (t).For very small feedback methane-temperature coefficients k we observe a linear growth methane concentration (see Fig. 2).For slightly larger coefficients we see a parabolic growth (see Fig. 3) and for larger ones one can see a true "Armageddon" (see Fig. 4).In the beginning, the contribution of H is 0.3 × 10 −3 part of βX , however, within hundred and thousand years this effect can be important.

Conclusion and discussion
Thawing permafrost and the resulting decomposition of previously frozen organic carbon to methane is one of the most significant potential feedback from terrestrial ecosystems to the atmosphere in a changing climate.Methane emissions from tundra permafrost lakes can produce a significant positive feedback of climate change.In this paper a new approach to modeling of methane emission from permafrost lakes is proposed.Usually the models of methane emission are based a number of differential equations describing heat and mass transfer in permafrost-climate system.These equations are complemented by number of empirical relations characterizing of a particular terrestrial ecosystems condition.While the approach is very reliable.It is difficult to be applied to complicated model for prediction of methane flux into the atmosphere on a global scale.On the other hand, we need to evaluate the positive feedback in

ESDD Introduction
Full the permafrost-climate system as a whole.We show how one can use macroscopic approaches in order to evaluate of this positive feedback.Indeed, in previous models it is required to use a number of factors to describe the microscopic processes (for example, diffusion of methane).These factors can only be set by in situ and only for each particular ecosystem.For realization of new macroscopic approach we can use meso-scale data such as an averaged lake radius.It gives us a new way: dynamics of permafrost lake can be calculated using remote sensing data; it does not require in situ measurements.In addition, the base on the proposed method we can estimate the coefficient of positive feedback of the permafrost-climate system and to predict possible methane flux into the atmosphere from thawing permafrost lakes.This paper is described only the test's estimate, which we can receive in terms of observational data and mathematical assumptions for phase transition theory.As a result, this research demonstrates that the methane emission can be occurring gradually as well as catastrophically.
Thus, incorporate the rigorous mathematical approach to different climate models (where the permafrost component is studied), we can describe the positive feedback for the permafrost-climate system in detail to come nearer to a theory of "Arctic Armageddon" (Kerr, 2010).Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | This equation can be improved and transformed to a more realistic model.It is clear that this growth is a complicated stochastic process.To take into account stochastic effects, we derive a Fokker-Planck equation corresponding to the deterministic model.This equation has, as a solution, the Pareto distribution.It is consistent with some Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

For
k = 0.8 we have simulated the Pareto distribution for lakes in Central Yakutia.It is known that this region has about 100 000 lakes with typical sizes R = 200 − 300 m, up to R = 1000 − 2000 m, that occupied the area about S 0 = 5 the standard algorithm of the Pareto law generation, give R min ≈ 80 m and A min ≈ 2000 m 2 , and the maximal lake radius 1.6 × 10 3 m, we see therefore that the Pareto law is, roughly, consistent with observed lake parameters.To estimate δ, we consider the area increment rate (S 1 − S 0 )/δT that gives δ ≈ 0.2 m per month.Notice that the termocarst lakes are a result of longtime evolution, thus, if Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Sudakov et al.
Finally we can assume that γ ∈ [0.2 × 10 resolve two elementary differential equations.The first is as follows dX dt = βX , (30) where the coefficient β = ln(1.01)defines the averaged observed methane growth.Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .Fig. 2 .
Fig. 1.Schematic of lake growth model.Blue curve: initial positions of thawing fronts; red curve: final positions of thawing fronts, which form a shall lake of bowl form; horizontal axe -horizontal radius; vertical axe: distance from surface to bottom.