## Abstract

We propose a two-dimensional model of a valley glacier in order to reconsider the question of whether thermal runaway could be a viable mechanism for the onset of creep instability in surging glaciers. We do this by providing an approximate solution for the temperature field based on the idea that shear is concentrated at the glacier bed. With this assumption, we show that a closed-form evolution equation for the glacier profile exists. While this is well known for isoviscous flows, it has not been previously derived for variable viscosity flows. During the process of deriving this equation, we show that thermal runaway does not occur. We provide numerical solutions of the model, and are led to infer that enhanced basal heating owing to refreezing of surface meltwater is an essential constituent in raising the bed temperature to the melting point.

## 1. Introduction

The explanation of the mechanism of why some glaciers surge has been a primary concern for glaciologists for many years. Meier & Post (1969) established many key properties of surging glaciers, and these have been variously explained by a variety of theories since. A glacier surge occurs when the glacier starts to slide rapidly at its base, and following the pioneering studies on Variegated Glacier in the 1980s (Kamb *et al*. 1985), it is now generally accepted that this rapid sliding is associated with subglacial water reaching elevated pressure. The same idea was indicated earlier by Lliboutry (1964, 1968) in his theoretical studies of glacier sliding.

One possible way in which high water pressures could be achieved is through the enhanced production of basal water, which in turn could be due to enhanced strain heating. This notion led to the suggestion that creep instability (Clarke *et al*. 1977) might be the cause of glacier surging, and in particular that the temperature in a glacier might undergo thermal runaway. Thermal runaway in shear flows of viscous fluids was studied by Gruntfest (1963) and Joseph (1964, 1965), and it was suggested as a possible causative surge mechanism by Clarke *et al*. (1977) and Yuen & Schubert (1979). These studies showed that multiple steady states were possible in glacier flow models of a prescribed depth, but their interpretation was criticized by Fowler (1980), on the basis that in reality one should prescribe the ice flux (via the net accumulation) and not the depth. Subsequently, Fowler & Larson (1980*a*,*b*) showed that, at least when thermal advection was negligible, a non-isothermal model for glacier flow indeed had a unique solution, and moreover this was linearly stable, suggesting that surge-like oscillations were unlikely in that case.

The neglect of thermal advection was unlikely to be critical to those studies because one would expect advection to have a stabilizing tendency. However, the possibility for thermal runaway could not be ruled out under conditions where the ice surface responds more slowly than the temperature field. In this case, thermal runaway could occur before the ice surface has time to respond.

The situation concerning thermal runaway is thus not entirely conclusive, and although it is not of active concern in valley glacier studies, the same issue arises in the study of ice streams, where the same precepts may be important (Payne & Dongelmans 1997; Hindmarsh 2009). Whereas in simple one-dimensional combustion studies one can prove the existence of thermal runaway (e.g. Fujita 1969), such precision is undoubtedly much harder in a glacier flow model involving a free boundary together with a nonlinear advection term. Our approach therefore is to seek an approximate theoretical framework that is both apparently conducive to the occurrence of thermal runaway, and susceptible to sufficient analysis that, within this framework, the issue can be satisfactorily resolved. While this falls short of an absolute proof one way or the other, it lends strong support to the more general validity of the conclusion, which is that runaway is unlikely to occur.

In this paper, we examine the possibility of thermal runaway by means of an approximation method based on the idea that shear is concentrated near the base of an ice flow (Nye 1959). Formally we assume that a parameter *γ* is large, indicating strong temperature dependence of the flow law. This allows us to solve the temperature equation analytically, and we are thus able to deduce an explicit evolution equation for the ice surface. Using this approach, we will show that thermal runaway is still not possible, and remains an unlikely cause of glacier surges.

## 2. Mathematical model for a valley glacier

The central focus of our study is the temperature equation, and we begin with this. As shown in figure 1, the downslope and transverse coordinates are *x* and *z*, with corresponding velocity components *u* and *w*. The absolute temperature is denoted by *T*, and satisfies the equation
2.1
where already we use the shallow ice approximation to neglect the longitudinal heat conduction term, on the basis that the aspect ratio of the flow (length divided by depth) is large. In this equation, *ρ* is the ice density, *c*_{p} is the specific heat, *k* is the thermal conductivity, *τ*_{ij} is the stress tensor and is the strain rate tensor. Again, the shallow ice approximation indicates that the only significant strain rate and stress are the shearing terms
2.2
where the second strain invariant *τ* satisfies
2.3
Glen’s flow law1
2.4
then implies the approximation
2.5

The shallow ice approximation for the momentum equations yields the expression
2.6
where *z*=*s* denotes the ice surface (we also take the bed to be at *z*=*b*), *χ* is the mean bed inclination angle and *g* is the acceleration due to gravity. Finally, the surface evolution equation follows from the consideration of conservation of mass:
2.7

We scale these equations by writing
2.8
where *T*_{m} is the melting temperature. The notation *x*∼*l* (and so on for the other variables) is the common abbreviation for the process of writing *x*=*l**x**, and then the dimensionless equations are written in terms of the starred variables; the stars are omitted for convenience. We could have equally written *T*−*T*_{m}∼Δ*T*, but it is useful to distinguish the scaled temperature with the different symbol *θ* as it has a different origin. We suppose that the rate factor *A* is thermally activated, thus
2.9
where
2.10
and we take Δ*T* to be the range of the prescribed surface temperature below the freezing point. Assuming Δ*T*≪*T*_{m}, we can write equation (2.9) in the approximate form
2.11
where
2.12

The resultant dimensionless equations are then
2.13
where the parameters are defined by
2.14
and we have used balance of mass conservation to choose *d* via
2.15

The boundary conditions are taken to be the prescription of a constant surface temperature −Δ*T* (Celsius), thus the dimensionless surface temperature is
2.16
The condition on the velocity is that of a basal sliding velocity *u*_{b}, at least if the basal temperature is at the melting point; however, if the basal temperature is below freezing, then we prescribe an effective geothermal heat flux *G**, and the basal velocity is then zero. The effective geothermal heat flux is the sum of the actual geothermal heat flux *G* and a term that represents the latent heat release by refreezing of surface (*not basal*) meltwater that finds its way to the base of the glacier through moulins and crevasses.2 If the quantity of ice melted, which heats the base in this way, is denoted *V* , measured as surface loss in elevation per year, then we have
2.17
where *L* is the latent heat. It turns out that this heat source dwarfs the geothermal heat flux in glaciers, and is instrumental in causing them to become warm at the base.

In dimensionless terms, these boundary conditions take the alternate forms, all applied at the base *z*=*b*,
2.18
where
2.19
is the basal stress, *u*_{b}(*τ*_{b}) is the fully temperate sliding law, and a function of the basal stress, and the dimensionless geothermal heat flux is
2.20

The four alternative conditions were provided by Fowler & Larson (1980*a*). The first three are referred to by them as cold, subtemperate and temperate, respectively, and we may term the last condition as warm; it refers to the situation where the ice above the bed is temperate and contains moisture. In this case, an enthalpy variable generalizes the temperature, and the dimensionless temperature *θ* is generalized to
2.21
where
2.22
and *ρ*_{w} is the density of water, *L* is the latent heat, *w* is the volumetric water fraction of temperate ice; St is the Stefan number. In the case that warm basal ice occurs, the energy equation needs to be modified to allow for the different physical process (Darcy flow) that allows the transport of enthalpy. This has been discussed by Fowler (2001), for example, but is not pursued here.

To derive estimates for the parameter values, we note that *u*_{0}, *τ*_{0} and *d* are not independent, but that from their definitions in equations (2.8) and (2.15), we have
2.23
We let
2.24
denote the elevation drop of the glacier, and note that then
2.25
The values of the parameters are given in table 1, and using these values we derive the estimates shown in table 2. We distinguish between two sets of values, recognizing that there is no one choice of values that can fit all locations. The primary difference lies in whether we suppose surface ablation occurs. In the normal case where this occurs (in subpolar and temperate climates), some of the surface meltwater will find its way to the bed, and if the base is frozen, some of this water is likely to refreeze, causing a release of latent heat. It may seem odd to allow such basal refreezing to occur at temperatures below the melting point: surely refreezing implies that the basal temperature is *at* the melting point. The point is resolved by realizing that one makes a distinction between the instantaneous temperature of an ice–water interface and the average temperature of the combined system. In a freezer, one can import trays of hot water and remove them when frozen, replacing them with more hot water. While the water is freezing, its interfacial temperature will be at the melting point, but the air temperature will be significantly lower. If the average temperature of the freezer and its contents is computed, it will remain below the melting point, and the net effect of the water input is to provide an effective source of heat.

The latent heat thus produced enhances the geothermal heat flux in a significant way. Where this does not occur (in polar, continental climate such as in Antarctica), there may be no melting at all, associated with a colder surface temperature (hence larger Δ*T*), and decreased precipitation (hence smaller *a*_{0}). These two climatic examples act as particular stereotypes in our investigations.

## 3. Approximation and reduction

None of the parameters in table 2 are particularly small or large, but a number of them conspire to suggest a useful means of approximation. The rate factor decreases (if *γ*=2.5) by a factor of 12 at the ice surface, and this decrease is enhanced by the stress dependence of the flow law. Half way to the surface, the creep rate is 3.5 times lower, but the consequent strain rate is then 28 times lower, because of the stress dependence. This effect is further enhanced by the relatively small value of *β*, which tends to concentrate the change of temperature (and thus also the shear) in a basal thermal boundary layer. Thus, the shear is concentrated at the glacier bed, and we can formally base an approximation scheme on this observation by considering the limit in which *γ*≫1, noting that the limits *β*→0 and enhance the quality of the approximation.

We write *γ**θ*=*ϕ*, so that the temperature equation takes the form
3.1
From this we see that *ϕ* will relax to equilibrium (if *s* is stationary) on a time scale of O(1/*α**γ*), and this steady state will be linear outside a boundary layer of thickness , because the exponential term becomes transcendentally small away from the boundary. Note that the advective term becomes small in this approximation, so that the quasi-steady state for *θ* is the solution of
3.2
together with the boundary conditions from equations (2.16) and (2.18).3 Noting that the viscous heating term is *α**τ**u*_{z}, we can integrate this equation once to find
3.3
and using equation (2.13)_{4}, we find that the evolution equation for *s* takes the form
3.4

We can progress further by noting that, as the exponential term in equation (3.2) is negligible away from the base, we can replace *τ* in that equation by *τ*_{b}, and the quasi-steady solution for *θ* can then be written explicitly as
3.5
where
3.6
and *B* must be positive in order that *θ* be real. The temperature gradient is calculated from equation (3.5), and is
3.7

To keep things simple, we now suppose that the bed of the glacier *b*=0, i.e. it is a flat inclined plane (figure 1), and we suppose also that the sliding velocity *u*_{b} is negligible. Thus, the subtemperate part of the base shrinks to a point, and the temperate and warm parts of the base satisfy the same condition. Thus, we suppose that *u*=0 at *z*=0 always, and
3.8
also at *z*=0.

The condition *θ*=−1 at *z*=*s* implies
3.9

### (a) Cold-based ice

It is well known (cf. Joseph 1964, 1965) that there are generally two solutions of equation (3.2) for *λ* less than some critical value *λ*_{c}; in the present case, these are distinguished as being a lower, cool branch, and an upper, warm branch.

For cold-based ice, we have in addition to equation (3.9) the equation
3.10
which shows that *C*>0. As we take *γ*≫1, equation (3.10) implies that *B*≪1, and therefore equation (3.9) is approximately
3.11
If we first suppose that *C*≤O(1), then
3.12
3.13
this solution therefore corresponds to the warm branch.

Alternatively, if *C*≫1, then
3.14
and in this case
3.15
and this is negative if *Γ**s*<1, which is also the criterion that *C*≫1. Thus, this cool lower branch is determined by equation (3.14), and thermal runaway occurs if *Γ**s*>1, except that in any case we then switch to the temperate-based case.

To compute , we use equation (3.7), together with the approximation
3.16
for large *ξ*. We then find that
3.17
and hence we find in this case that *s* satisfies
3.18

Because, particularly in the subpolar case, *γ* is not in fact that large, one must be careful that the approximations in equation (3.14) are reasonably accurate. Figure 2 shows two examples of the exact solutions for *B* as a function of *s* computed from the solution of equations (3.9) and (3.10), together with the approximation . In each curve, the upper branch is the stable cool branch, and one can see that the approximate value actually gives a good estimate, even though is not so very small. This gives us some confidence in the use of equation (3.14).

### (b) Temperate-based ice

For the case that *θ*=0 at *z*=0, we have instead of equation (3.10)
3.19
From equations (3.7) and (3.5), we see that if *C*>0, then *θ*_{z}|_{0}<0, and *θ*<0 in *z*>0, whereas if *C*<0, the opposite is true, and the basal ice is warm (actually, superheated). In addition, the cool branch is stable and the warm branch is unstable.

We can see from equation (3.19) that either *B*=O(1) or *B*≪1. If *B*=O(1), then equation (3.9) implies (as *γ*≫1) that |*C*|≫1, and hence in fact *B*≪1 is the only possibility. The two possible solutions then correspond to ±*C*≫1: for the stable, cool branch, *C*≫1, while for the unstable, warm branch, −*C*≫1. From equation (3.19), we have
3.20
and therefore, from equation (3.9), we obtain the approximations
3.21
Restricting our attention to the cool, stable case *C*>0, we can now approximate as before, using the definitions in equations (3.7) and (3.21), and the approximation (3.16). The result is
3.22
Finally, using the definition of *λ* in equation (3.6), we find that for temperate basal conditions, the mass conservation equation (3.4) takes the form
3.23

The temperate-based evolution equation applies while , that is,
3.24
The warm and cool branches must come together when *C*=O(1), and thus *B*=O(1). Thermal runaway will thus occur if *B*=O(1). Formally this does not happen, as *B*≪1. If we suppose the approximate result in equation (3.21) can be extended, then it suggests that runaway will occur if for some value of *B*=O(1), or
3.25
but just as before, we can associate the transition from cool to warm branches with the passage of the basal heat flux, −*θ*_{z}|_{0}, through zero, and thus the onset of a region of basal moist ice; again, thermal runaway does not occur.

### (c) Summary

Both cold-and temperate-based glacier models can be represented by the combined form
3.26
where *K*(*s*) can be defined by
3.27
Note that *K* is a continuous function of *s*. The terms e^{−γ} are formally negligible, but are retained for the purpose of indicating the appropriate continuity at *Γ**s*=1.

In computing the temperature field, we would use the approximation given by equation (3.5), where the values of *B* and *C* should be as described above. With the approximate formulae given by equations (3.14) and (3.21)_{1}, the temperature field will be continuous across the cold–temperate transition at *Γ**s*=1, but the basal temperature will not be exactly zero at the switch point because equation (3.19) is only approximately satisfied there. Because it represents the actual melting point we want to ensure that the transition occurs at *θ*|_{0}=0, and the simplest way to do this, if we want to retain explicit approximations for *B* and *C*, rather than solving for them exactly at each value of *x*, is to make the same approximation in equation (3.5) as we do in deriving the approximate values of *B* and *C*.

As this involves the supposition that *B* is small and *C* is large, we make the same approximation in equation (3.5), and after some algebra, this leads to the linear approximation
3.28

## 4. Numerical results

In deriving the effective non-isothermal surface evolution profile equation (3.26), we have tacitly assumed *s*_{x}<1/*μ*, which physically implies the glacier surface is inclined downhill, as we expect. More generally, the derivation of the equation leads to the form
4.1

### (a) Boundary conditions

The issue of boundary conditions for a valley glacier is an interesting one. The second-order equation (4.1) requires two boundary conditions, upstream and downstream. That at the snout of the glacier, assuming it terminates on land, is simply
4.2
where *x*_{f} is the position of the glacier front. The value of *x*_{f} is unknown (it is a free boundary), but the condition (4.2) is sufficient to determine it, as the equation is degenerate there (the diffusion coefficient is zero).

At the upstream end of the glacier *x*=0, the position is less clear. We might suggest
4.3
but the flux
4.4
is then apparently zero unless *s*_{x} is infinite. Indeed, the condition *s*=0 then requires *q*∼*a**x*>0 near *x*=0, and this requires *s*_{x} to be infinite and *negative*. Thus, the only way one can maintain *s*=0 at the head is to have a finite flux, and (since *s* must be positive in *x*>0) this must be a negative flux; that is to say, the glacier actually grows into an ice sheet with an ice divide downstream of the head.

Physically, we need to represent the bergschrund (that is to say, the head) of the glacier. One way to do this is to allow a variable basal slope, thus replacing equation (4.1) by
4.5
where *Θ* is the scaled basal slope. We can represent a bergschrund by specifying
4.6
where 0<*ν*<1 in order that the mountain height be finite. Then, we can have positive downstream flux with *q*∼*a**x* near *x*=0 if *s*∼*x*^{σ}, where
4.7

While this provides an algebraically convenient resolution to the problem, and allows us to prescribe equation (4.3), it introduces unnecessary numerical complication. Thus, in practice, we have applied a flux boundary condition
4.8
With a small positive value for *q*_{0}, then *s*_{x}<1/*μ*, and *s* is finite at *x*=0.

### (b) Results

In this section, we present numerical results corresponding to three particular cases: cold, temperate and polythermal. The results have been obtained by using the numerical methods briefly described in the appendix.

The departure point for the numerical techniques is a fixed domain formulation for the moving boundary problem. Particularly, we select the spatial domain *Ω*=[0,*M*], where *M* is sufficiently large that *Ω* always contains the glacier profile. In practice, we take *M*=3. Next, we fix the boundary conditions at both ends of the domain. As intimated above, these are
4.9
where *q* is given by equation (4.4) (and *K* by equation (3.27)). The initial condition is taken as
4.10
In all runs, we have chosen the accumulation rate function to be
4.11
and the Glen exponent *n*=3.

For the numerical solution, a spatial mesh stepsize Δ*x*=10^{−4} and a time step of Δ*t*=10^{−6} have been used. Moreover, for the fixed-point iteration to deal with the nonlinear diffusive term, a stopping test based on the error between two iterations being less than 10^{−6} has been used, while for the duality method 10^{−9} has been chosen (see the appendix for further details on numerical methods).

### (c) Example 1. Cold-based ice

The first example corresponds to the polar parameters that may be appropriate for the case where the effective basal heating is very small, and we expect to find a cold-based glacier. The parameters used are
4.12
Figure 3 shows the computed profile evolution from *t*=1 to *t*=4 (when the profile has reached the steady state), after numerically solving equation (3.26). Figure 4 shows the computed temperature for *t*=4 in the ice region, obtained from equation (3.28).

### (d) Example 2. Temperate-based ice

The second example corresponds to the subpolar parameters for which the effective basal melting is large, and we expect a temperate-based glacier.4 The parameters used are
4.13
In figure 5, we show the computed profile evolution from *t*=1 to *t*=4 (when it has again reached steady state), after numerically solving equation (3.26). Figure 6 shows the computed temperature for *t*=4 in the ice region, obtained from equation (3.28).

### (e) Example 3. Polythermal basal ice

The final example corresponds to a case where basal heat flux is smaller than for the subpolar parameters. We expect that as *Γ* is reduced, the base will become increasingly cold at the head and snout, and that eventually the glacier will become cold based. Between the temperate-based and cold-based states, we find a mixed cold–temperate, or polythermal basal régime at the base of the glacier. The parameters used are
4.14
In figure 7, we show the computed profile evolution from *t*=1 to *t*=4 (again at steady state), after numerically solving equation (3.26). Figure 8 shows the computed temperature for *t*=4 in the ice region, obtained from equation (3.28).

## 5. Conclusions

Our principal aim in this paper was to consider the issue of whether thermal runaway could properly happen in glacier flow, in a more comprehensive way than has yet been done. In doing so, there are two problems to consider. The first is to consider whether runaway can occur in the different possible basal thermal régimes that can exist at the base of a glacier, and the second is to consider whether, if runaway does occur, it leads, in practice, to an unlimited acceleration within the same basal thermal régime. Earlier work has shown that thermal runaway in shear flow models with fixed depth can occur, and this therefore remains a possibility if the temperature equation can respond more rapidly than the ice surface. Formally, this is the case for strongly temperature-dependent rheology, here taken as the limit *γ*≫1.

In indulging this limit, we therefore provide a glacier flow model that would seem most likely to produce runaway, if it can occur at all. However, we have found that the progressive adjustment of thermal boundary conditions from cold to subtemperate to temperate never in practice allows runaway to occur, because one simply switches basal boundary conditions at the point where runaway could occur. The possibility of runaway when the basal *ice* (not boundary condition) becomes temperate, and internal moisture is generated, remains a possibility, but realistic modelling of temperate ice, allowing for moisture drainage, has not been done.

In building our approximate model to address this problem, we have constructed an effective non-isothermal ice surface evolution equation (4.1) that represents in a realistic way the evolution of the ice surface when the flow law is temperature dependent. The derivation of this equation suggests that temperature variation has very little effect on the overall motion of the glacier, other than mildly adjusting the depth scale of the flow.

A particular revelation of this exercise has been the unexpected inadequacy of geothermal heat. As measured by the parameter *Γ*, we have shown that a glacier remains cold based until its depth *s*>1/*Γ*. If *Γ* is computed using normal values of geothermal heat flux, we find that it is so small that basal ice will never reach the melting point. In reality, entirely cold-based glaciers are something of a rarity, and we consider this to be due to the overwhelming importance of latent heat release by buried surface meltwater and rainwater. In the absence of such enhanced basal heating, glaciers would remain frozen at their base.

## Acknowledgements

A.C.F. acknowledges the support of the Mathematics Applications Consortium for Science and Industry (www.macsi.ul.ie) funded by the Science Foundation Ireland mathematics initiative grant 06/MI/005. C.V. and R.T. acknowledge the support of research projects of MCII (MTM2007–67596–CO2–01) and Xunta de Galicia (PGIDIT06PX105230PR).

## Appendix: Numerical methods

In this section, we briefly describe the set of numerical methods that have been used for the three models. We write equation (4.1) in the form
A1
where
A2
and *K* is given by equation (3.27).

Equation (A1) is valid for the function *s* when *s*(*x*,*t*)>0, the set of whose points is an additional unknown. Thus, following previous papers (e.g. Calvo *et al.* 2000), we use a fixed domain technique.

For this purpose, let ; in our computations we use *M*=3. If we extend the function *s* to be zero in the ice-free region, then the extended function (still noted by *s*) satisfies the equations
A3
Notice that in the ice covered region, the fact that *s*(*x*,*t*)>0 together with equation (A3)_{3} implies equation (A1).

Next, in order to solve numerically the complementarity problem (A3)_{1−4}, an implicit finite-difference scheme for time discretization is applied together with an iterative fixed-point technique for the resulting nonlinear diffusive term.

So, for each , we initialize *s*^{m+1,0} and the following problem has to be solved at step *k*+1:

find *s*^{m+1,k+1} such that
A4
where
A5
In order to pose the variational formulation and the spatial discretization, we point out that the relation between variational inequalities, complementarity problems and obstacle problems can be reviewed in Elliott & Ockendon (1981), for example.

Next, in order to discretize in space the variational inequality associated with equation (A4), a piecewise linear Lagrange finite-element space is used. Thus, for a given positive mesh stepsize Δ*x*, a uniform finite-element mesh *τ*_{Δx} is built for the domain *Ω* with nodes *x*_{i}=(*i*−1)Δ*x*, *i*=1,…,*N*+1. The following classic spaces and sets are introduced:
A6
where *E* denotes a standard finite element. In this way, the discretized problem can be written as follows:

find such that
A7
In order to solve the discretized nonlinear problem (A7), several minimization algorithms applied to variational formulations of obstacle problems can be used (Glowinski 1984). In this work, a duality method proposed in Bermúdez & Moreno (1981) to solve variational inequalities has been chosen. This method has been successfully applied previously in obstacle problems arising in glaciology (e.g. Calvo *et al*. (2000) for further technical details).

## Footnotes

↵1 Usually taken with

*n*=3, which we also do in this paper.↵2 This condition does not imply that the basal temperature is at the melting point: think of the temperature regulation of a freezer where one regularly adds ice cube trays containing water.

↵3 There is a slight inaccuracy here. The advective term is certainly small in the boundary layer, but both the exponential and conductive terms become small in the bulk flow, so that in fact the outer problem for

*ϕ*is the approximate equation**u**⋅**∇***ϕ*≈0, if we suppose*β*≪1, the solution to which is*ϕ*=*ϕ*(*ψ*), where*ψ*is a stream function for the flow. With the constant surface temperature condition (2.16), this makes no difference at all, and in fact even with a varying surface temperature, the boundary layer solution is formally the same, provided we take the dimensionless surface temperature*θ*=−1 to be that at the head*x*=0.↵4 Over most of the domain. Near the snout of the glacier, the base must become cold according to the model.

- Received July 26, 2009.
- Accepted September 25, 2009.

- © 2009 The Royal Society