Mathematics of Planet Earth at Banff

7 September, 2011

A while back, I mentioned that 2013 will be a special year for programs on the Mathematics of Planet Earth. I also mentioned that the Centre de Recerca Matematica in Barcelona is inviting mathematicians to organize conferences and workshops on this theme.

They’re also inviting mathematicians to organize workshops on this theme at the Banff International Research Station for Mathematical Innovation and Discovery, or BIRS. This is a famous and beautiful research center in the Canadian Rockies.

The deadline is coming up on September 30th, and I want to apply. If you’d like to join me, please drop me a note, either here on this blog or by email!

I’m open to all sorts of ideas, and I’d love help from biologists or climate scientists. If you don’t give me a better idea, I’ll probably do an application on network theory. It might look a bit like this:

Diagrammatic languages for describing complex networks made of interacting parts are used throughout ecology, biology, climate science, engineering, and many other fields. Examples include Systems Biology Graphical Notation, Petri nets in computer science, stochastic Petri nets and chemical reaction networks in chemistry and biochemistry, bond graphs in electrical, chemical and mechanical engineering, Bayesian networks in probabilistic reasoning, box models in climate science, and Harold Odum’s Energy Systems Language for systems ecology. Often these diagrammatic languages are invented by practitioners in a given field without reference to previous work in other fields. Recently mathematicians have set up the theoretical infrastructure needed to formalize, rigorously relate, and some cases unify these various languages. Doing this will help interdisciplinary work of the sort that is becoming important in theoretical ecology, climate science and ‘the mathematics of planet Earth’. The goal of this workshop is to bring together experts on various diagrammatic languages and mathematicians who study the general theory of diagrammatic reasoning.

If you’d be interested in coming to a workshop on this subject, let me know. Banff provides accommodation, full board, and research facilities—but not, I believe, travel funds! So, “interested in coming” means “interested enough to pay for your own flight”.

Banff does “full workshops” with 42 people for 5 days, and “half workshops” with 20 people for 5 days. Part of why I’m asking you to express your interest is to gauge which seems more appropriate.

Here’s what they say:

With a growing global population competing for the same global resources, an increased frequency and intensity of dramatic climatic events, and evidence pointing to more long-term patterns of general climate change, the pressure to comprehend nature and its trends is greater than ever. Leaders in politics, sociology and economics have begun to seriously take note of issues which before were confined to the natural sciences alone, and mathematical modeling is at the heart of much of the research undertaken. The year 2013 has thus been earmarked by mathematical sciences institutes around the world as a time for a special emphasis on the study of the “Mathematics of Planet Earth” (MPE 13). This theme is to be interpreted as broadly as possible, in the aim of creating new partnerships with related disciplines and casting new light on the many ways in which the mathematical sciences can help to comprehend and tackle some of the world’s most pressing problems.

The Banff International Research Station (BIRS) is a full partner in this important initiative, as the goals of MPE 13 are completely in line with the station’s commitment to pursuing excellence in a broad range of mathematical sciences and applications. BIRS has already planned to host five workshops in 2012 which deal with the themes of MPE 13:

“Emergent Behavior in Multi-particle Systems with Non-local Interactions” (January 22-27).

“Frontiers in the Detection and Attribution of Climate Change” (May 29–June 1).

“Tissue Growth and Morphogenesis: from Genetics to Mechanics and Back” (July 22-27).

“Model Reduction in Continuum Thermodynamics: Modeling, Analysis and Computation” (September 16-21).

“Thin Liquid Films and Fluid Interfaces: Models, Experiments and Applications” (December 9-14).

BIRS also invites interested applicants to use the opportunities of its 2013 program and submit proposals in line of the MPE 2013 theme, in conjunction with BIRS’ regular format for programming. Proposals should be made using the BIRS online submission process.


US Weather Disasters in 2011

6 September, 2011

The US Federal Emergency Management Agency (FEMA) is running out of money!

So far this year, ten weather disasters have each caused over a billion dollars of damage in the United States. This beats the record set in 2008, when there were nine. FEMA now has less than a billion dollars in its fund:

• Brian Naylor, Costs Of Irene Add Up As FEMA Runs Out Of Cash, Morning Edition, National Public Radio, 30 August 2011.

Let’s review these disasters:

10) Hurricane Irene, August 27-28: A large and powerful Atlantic hurricane that left extensive flood and wind damage along its path through the Caribbean, the east coast of the US and as far north as Atlantic Canada. Early estimates say Irene caused $7 billion in damages in the US.

9) Upper Midwest flooding, summer: An above-average snowpack across the northern Rocky Mountains, combined with rainstorms, caused the Missouri and Souris rivers to swell beyond their banks across the Upper Midwest. An estimated 11,000 people were forced to evacuate Minot, N.D. Numerous levees were breached along the Missouri River, flooding thousands of acres of farmland. Over $2 billion in damages.

8) Mississippi River flooding, spring-summer: Persistent rainfall (nearly triple the normal amount in the Ohio Valley), combined with melting snowpack, caused historical flooding along the Mississippi River and its tributaries. At least two people died. $2 to $4 billion in damages.

7) Southern Plains/Southwest drought, heat wave and wildfires, spring and summer: Drought, heat waves, and wildfires hit Texas, Oklahoma, New Mexico, Arizona, southern Kansas, western Arkansas and Louisiana this year. Wildfire fighting costs for the region are about $1 million per day, with over 2,000 homes and structures lost by mid-August. Over $5 billion in damages so far.

6) Midwest/Southeast tornadoes, May 22-27: Central and southern states saw approximately 180 twisters and 177 deaths within a week. A tornado in Joplin, Mo., caused at least 141 deaths—the deadliest single tornado to strike the United States since modern record keeping began in 1950. Over $7 billion in damages.

5) Southeast/Ohio Valley/Midwest tornadoes, April 25-30: This outbreak of tornadoes over central and southern states led to 327 deaths. Of those fatalities, 240 occurred in Alabama. The deadliest of the estimated 305 tornadoes in the outbreak was an EF-5 that hit northern Alabama, killing 78 people. Several big cities were directly affected by strong tornadoes, including Tuscaloosa, Birmingham and Huntsville in Alabama, and Chattanooga in Tennessee. Over $9 billion in damages.

4) Midwest/Southeast tornadoes, April 14-16: An outbreak over central and southern states produced an estimated 160 tornadoes. Thirty-eight people died, 22 of them in North Carolina. Over $2 billion in damages.

3) Southeast/Midwest tornadoes, April 8-11: An outbreak of tornadoes over central and southern states saw an estimated 59 tornadoes. Over $2.2 billion in damages.

2) Midwest/Southeast tornadoes, April 4-5: An outbreak of tornadoes over central and southern states saw an estimated 46 tornadoes. Nine people died. Over $2.3 billion in damages.

1) Blizzard, Jan 29-Feb 3: A large winter storm hit many central, eastern and northeastern states. 36 people died. Over $2 billion in damages.

I got most of this information from this article, which was written before Irene pushed 2011 into the lead:

• Brett Israel, 2011 ties for most billion-dollar weather disasters, Our Amazing Planet, 18 August 2011.

We can expect more weather disasters as global warming proceeds. The National Academy of Sciences says:

• Increases of precipitation at high latitudes and drying of the already semi-arid regions are projected with increasing global warming, with seasonal changes in several regions expected to be about 5-10% per degree of warming. However, patterns of precipitation show much larger variability across models than patterns of temperature.

• Large increases in the area burned by wildfire are expected in parts of Australia, western Canada, Eurasia and the United States.

• Extreme precipitation events—that is, days with the top 15% of rainfall—are expected to increase by 3-10% per degree of warming.

• In many regions the amount of flow in streams and rivers is expected to change by 5-15% per degree of warming, with decreases in some areas and increases in others.

• The total number of tropical cyclones should decrease slightly or remain unchanged. Their wind speed is expected to increase by 1-4% per degree of warming.

Some people worry about sea level rise, but I think the bite from weather disasters and ensuing crop failures will hurt much more, much sooner.

Since it doesn’t look like politicians will do enough to cut carbon emissions, insurance companies are moving to act on their own—not to prevent weather disasters, but to minimize their effect:

Swiss Re’s global headquarters face Lake Zurich, overlooking a small yacht harbor. Bresch and a colleague, Andreas Schraft, sometimes walk the 20 minutes to the train station together after work, past more yachts, an arboretum, and a series of bridges. In September 2005, probably on one of these walks, the two began to discuss what they now call “Faktor K,” for “Kultur”: the culture factor. Losses from Hurricanes Katrina, Rita, and Wilma had been much higher than expected in ways the existing windstorm models hadn’t predicted, and it wasn’t because they were far off on wind velocities.

The problem had to do more with how people on the Gulf Coast were assessing windstorm risk as a group. Mangrove swamps on the Louisiana coast had been cut down and used as fertilizer, stripping away a barrier that could have sapped the storm of some of its energy. Levees were underbuilt, not overbuilt. Reinsurers and modeling firms had focused on technology and the natural sciences; they were missing lessons from economists and social scientists. “We can’t just add another bell and whistle to the model,” says Bresch, “It’s about how societies tolerate risk.”

“We approach a lot of things as much as we can from the point of statistics and hard data,” says David Smith, head of model development for Eqecat, a natural hazards modeling firm. “It’s not the perfect expression.” The discrepancy between the loss his firm modeled for Katrina and the ultimate claims-based loss number for his clients was the largest Smith had seen. Like others in the industry, Eqecat had failed to anticipate the extent of levee failure. Construction quality in the Gulf states before Katrina was poorer than anticipated, and Eqecat was surprised by a surge in demand after the storm that inflated prices for labor and materials to rebuild. Smith recognizes that these are questions for sociologists and economists as well as engineers, and he consults with the softer sciences to get his models right. But his own market has its demands, too. “The more we can base the model on empirical data,” he says, “the more defendable it is.”

After their walk around the lake in 2005, Swiss Re’s Bresch and Schraft began meeting with social scientists and laying out two goals. First, they wanted to better understand the culture factor and, ultimately, the risks they were underwriting. Second, they wanted to use that understanding to help the insured prevent losses before they had to be paid for.

The business of insurers and reinsurers rests on balancing a risk between two extremes. If the risk isn’t probable enough, or the potential loss isn’t expensive enough, there’s no reason for anyone to buy insurance for it. If it’s too probable and the loss too expensive, the premium will be unaffordable. This is bad for both the insured and the insurer. So the insurance industry has an interest in what it calls “loss mitigation.” It encourages potential customers to keep their property from being destroyed in the first place. If Swiss Re is trying to affect the behavior of the property owners it underwrites, it’s sending a signal: Some behavior is so risky that it’s hard to price. Keep it up, and you’ll have no insurance and we’ll have no business. That’s bad for everyone.

To that end, Swiss Re has started speaking about climate risk, not climate change. That the climate is changing has been established in the eyes of the industry. “For a long time,” says Bresch, “people thought we only needed to do detailed modeling to truly understand in a specific region how the climate will change. … You can do that forever.” In many places, he says, climate change is only part of the story. The other part is economic development. In other words, we’re building in the wrong places in the wrong way, so wrong that what we build often isn’t even insurable. In an interview published by Swiss Re, Wolf Dombrowsky, of the Disaster Research Center at Kiel University in Germany, points out that it’s wrong to say that a natural disaster destroyed something; the destruction was not nature’s fault but our own.

In 1888 the city of Sundsvall in Sweden, built of wood, burned to the ground. A group of reinsurers, Swiss Re among them, let Sweden’s insurers know there was going to be a limit in the future on losses from wooden houses, and it was going to be low. Sweden began building with stone. Reinsurance is a product, but also a carrot in the negotiation between culture and reality; it lets societies know what habits are unsustainable.

More recently, the company has been working with McKinsey & Co., the European Commission, and several environmental groups to develop a methodology it calls the “economics of climate adaptation,” a way to encourage city planners to build in a way that will be insurable in the future. A study of the U.K. port of Hull looks at potential losses by 2030 under several different climate scenarios. Even under the most extreme, losses were expected to grow by $17 million due to climate change and by $23 million due to economic growth. How Hull builds in the next two decades matters more to it than the levels of carbon dioxide in the air. A similar study for Entergy (ETR), a New Orleans-based utility, concluded that adaptations on the Gulf Coast—such as tightening building codes, restoring wetlands and barrier islands, building levees around chemical plants, and requiring that new homes in high-risk areas be elevated—could almost completely offset the predicted cost of 100-year storms happening every 40 years.

I actually disagree somewhat with the statement “it’s wrong to say that a natural disaster destroyed something; the destruction was not nature’s fault but our own.” There’s some truth to this, but also some untruth. The question of “fault” or “blame” is a slippery one here, and there’s probably no way to completely settle it.

Is it the “fault” of people in Vermont that they weren’t fully prepared for a hurricane? After all, it’s rare—or at least it used to be rare—for hurricanes to make it that far north. The governor of Vermont, Peter Shumlin, recently said:

I find it extraordinary that so many political leaders won’t actually talk about the relationship between climate change, fossil fuels, our continuing irrational exuberance about burning fossil fuels, in light of these storm patterns that we’ve been experiencing.

We had storms this spring that flooded our downtowns and put us through many of the same exercises that we’re going through right now. We didn’t used to get weather patterns like this in Vermont. We didn’t get tropical storms. We didn’t get flash flooding.

We in the colder states are going to see the results of climate change first. Myself, Premier Charest up in Quebec, Governor Cuomo over in New York, we understand that the flooding and the extraordinary weather patterns that we’re seeing are a result of our burnings of fossil fuel. We’ve got to get off fossil fuels as quickly as we know how, to make this planet livable for our children and our grandchildren.

On the other hand, you could say that it is the fault of Vermonters, or at least humanity as a whole, for causing global warming in the first place.

But ultimately, pinning blame on someone or something is less important than figuring out how to solve the problems we face.


Melting Permafrost (Part 1)

1 September, 2011

Some people worry about rising sea levels due to global warming. But that will happen slowly. I worry about tipping points.

The word “tipping point” should remind you of pushing on a glass of water. If you push it a little, and then stop, it’ll right itself: no harm done. But if you push it past a certain point, it starts tipping over. Then it’s hard to stop.

So, we need to study possible tipping points in the Earth’s climate system. Here’s a list of them:

Tipping point, Azimuth Library.

Today I want to talk about one: melting permafrost. When melting permafrost in the Arctic starts releasing lots of carbon dioxide and methane—a vastly more potent greenhouse gas—the Earth will get even hotter. That, in turn, will melt even more permafrost. In theory, this feedback loop could tip the Earth over to a much hotter state. But how much should we worry about this?

Climate activist Joe Romm takes it very seriously:

• Joe Romm, NSIDC bombshell: Thawing permafrost feedback will turn Arctic from carbon sink to source in the 2020s, releasing 100 billion tons of carbon by 2100, Climate Progress, 17 February 2011.

If you click on just one link of mine today, let it be this! He writes in a clear, snappy way. But let me take you through some of the details in my own more pedestrian fashion.

For starters, the Arctic is melting. Here’s a graph of Arctic sea ice volume created by the Pan-Arctic Ice Ocean Modeling and Assimilation System—click to enlarge:

The blue line is the linear best fit, but you can see it’s been melting faster lately. Is this a glitch or a new trend? Time will tell.

2011 is considerably worse than 2007, the previous record-holder. Here you can clearly see the estimated total volume in thousands of cubic kilometers, and how it changes with the seasons:


As the Arctic melts, many things are changing. The fabled Northwest Passage is becoming a practical waterway, so battles are starting to heat up over who controls it. The U.S. and other nations see it as an international waterway. But Canada says they own it, and have the right to regulate and protect it:

• Jackie Northam, Arctic warming unlocking a fabled waterway, Morning Edition, National Public Radio, 15 August 2011.

But the 800-pound gorilla in the room is the melting permafrost. A lot of the Arctic is covered by permafrost, and it stores a lot of carbon, both as peat and as methane. After all, peat is rotten plant material, and rotting plants make methane. Recent work estimates that between 1400 and 1700 gigatonnes of carbon is stored in permafrost soils worldwide:

• C. Tarnocai, J. G. Canadell, E. A. G. Schuur, P. Kuhry, G. Mazhitova, and S. Zimov, Soil organic carbon pools in the northern circumpolar permafrost region, Global Biogeochemical Cycles 23 (2009), GB2023.

That’s more carbon than currently resides in all living things, and twice as much carbon as held by the atmosphere!

How much of this carbon will be released as the Arctic melts—and how fast? There’s a new paper about that:

• Kevin Schaefer, Tingjun Zhang, Lori Bruhwiler, Andrew Barrett, Amount and timing of permafrost carbon release in response to climate warming, Tellus B 63 (2011), 165–180.

It’s not free, but you can read Joe Romm’s summary. Here’s their estimate on how carbon will be released by melting permafrost:

So, they’re guessing that melting permafrost will release a gigatonne of carbon per year by the mid-2030s. Moreover, they say:

We predict that the PCF [permafrost carbon feedback] will change the Arctic from a carbon sink to a source after the mid-2020s and is strong enough to cancel 42-88% of the total global land sink. The thaw and decay of permafrost carbon is irreversible and accounting for the PCF will require larger reductions in fossil fuel emissions to reach a target atmospheric CO2 concentration.

One of the authors explains more details here:

“The amount of carbon released [by 2200] is equivalent to half the amount of carbon that has been released into the atmosphere since the dawn of the industrial age,” said NSIDC scientist Kevin Schaefer. “That is a lot of carbon.”

The carbon from permanently frozen ground known as permafrost “will make its impact, not only on the climate, but also on international strategies to reduce climate change Schaefer said. “If we want to hit a target carbon concentration, then we have to reduce fossil fuel emissions that much lower than previously calculated to account for this additional carbon from the permafrost,” Schaefer said. “Otherwise we will end up with a warmer Earth than we want.”

The carbon comes from plant material frozen in soil during the ice age of the Pleistocene: the icy soil trapped and preserved the biomass for thousands of years. Schaefer equates the mechanism to storing broccoli in the home freezer: “As long as it stays frozen, it stays stable for many years,” he said. “But you take it out of the freezer and it will thaw out and decay.”

Now, permafrost is thawing in a warming climate and “just like the broccoli” the biomass will thaw and decay, releasing carbon into the atmosphere like any other decomposing plant material, Schaefer said. To predict how much carbon will enter the atmosphere and when, Schaefer and coauthors modeled the thaw and decay of organic matter currently frozen in permafrost under potential future warming conditions as predicted by the Intergovernmental Panel on Climate Change.

They found that between 29-59 percent of the permafrost will disappear by 2200. That permafrost took tens of thousands of years to form, but will melt in less than 200, Schaefer said.

Sound alarmist? In fact, there are three unrealistically conservative assumptions built into this paper:

1) The authors assume the ‘moderate warming’ scenario called A1B, which has atmospheric concentrations of CO2 reaching 520 ppm by 2050 and stabilizing at 700 ppm in 2100. But so far we seem to be living out the A1F1 scenario, which reaches 1000 ppm by century’s end.

2) Their estimate of future temperatures neglects the effect of greenhouse gases released by melting permafrost.

3) They assume all carbon emitted by permafrost will be in the form of CO2, not methane.

Point 2) means that the whole question of a feedback loop is not explored in this paper. I understand why. To do that, you can’t use someone else’s climate model: you need to build your own! But it’s something we need to study. Do you know anyone who is? Joe Romm says:

Countless studies make clear that global warming will release vast quantities of greenhouse gases into the atmosphere this decade. Yet, no climate model currently incorporates the amplifying feedback from methane released by a defrosting tundra.

If we try to understand this feedback, point 3) becomes important. After all, while methane goes away faster than CO2, its greenhouse effect is much stronger while it lasts. For the first 20 years, methane has about 72 times the global warming potential of carbon dioxide. Over the first 100 years, it’s about 25 times as powerful.

Let’s think about that a minute. In 2008, we burnt about 8 gigatonnes of carbon. If Schaefer et al are right, we can expect 1 extra gigatonne of carbon to be released from Arctic permafrost by around 2035. If that’s almost all in the form of carbon dioxide, it makes our situation slightly worse. But if a lot of it is methane, which is—let’s roughly say—72 times as bad—then our situation will be dramatically worse.

But I don’t know how much of the carbon released will be in the form of methane. I also don’t know how much of the methane will turn into other organic compounds before it gets into the atmosphere. I’d really like to know!

I hope you learn more about this stuff and help me out. Here are a few good references available for free online, to get started:

• Edward A. G. Schuur et al, Vulnerability of permafrost carbon to climate change: implications for the global carbon cycle, Bioscience 58 (2008), 701-714.

• David M. Lawrence, Andrew G. Slater, Robert A. Tomas, Marika M. Holland and Clara Deser, Accelerated Arctic land warming and permafrost degradation during rapid sea ice loss, Geophysical Research Letters 35 (2008), L11506.

• Amanda Leigh Mascarelli, A sleeping giant?, Nature Reports Climate Change, 5 March 2009.

The last one discusses the rise in atmospheric methane that was observed in 2007:

It also discusses the dangers of methane being released from ice-methane crystals called methane clathrates at the bottom of the ocean—something I’m deliberately not talking about here, because it deserves its own big discussion. However, there are also clathrates in the permafrost. Here’s a picture by W. F. Kuhs, showing what methane clathrate looks like at the atomic scale:

The green guy in the middle is methane, trapped in a cage of water molecules. Click for more details.

If you know more good references, please tell me about them here or add them to:

Permafrost, Azimuth Library.


Hierarchical Organization and Biological Evolution (Part 1)

29 August, 2011

guest post by Cameron Smith

Consider these quotes:

My thesis has been that one path to the construction of a non-trivial theory of complex systems is by way of a theory of hierarchy. Empirically, a large proportion of the complex systems we observe in nature exhibit hierarchic structure. On theoretical grounds we could expect complex systems to be hierarchies in a world in which complexity had to evolve from simplicity.Herbert Simon, 1962

(Many of the concepts that) have dominated scientific thinking for three hundred years, are based upon the understanding that at smaller and smaller scales—both in space and in time—physical systems become simple, smooth and without detail. A more careful articulation of these ideas would note that the fine scale structure of planets, materials and atoms is not without detail. However, for many problems, such detail becomes irrelevant at the larger scale. Since the details (become) irrelevant (at such larger scales), formulating theories in a way that assumes that the detail does not exist yields the same results as (theories that do not make this assumption).Yaneer Bar-Yam

Thoughts like these lead me to believe that, as a whole, we humans need to reassess some of our approaches to understanding. I’m not opposed to reductionism, but I think it would be useful to try to characterize those situations that might require something more than an exclusively reductionist approach. One way to do that is to break down some barriers that we’ve constructed between disciplines. So I’m here on Azimuth trying to help out this process.

Indeed, Azimuth is just one of many endeavors people are beginning to work on that might just lead to the unification of humanity into a superorganism. Regardless of the external reality, a fear of climate change could have a unifying effect. And, if we humans are simply a set of constituents of the superorganism that is Earth’s biosphere, it appears we are its only candidate germ line. So, assuming we’d like our descendants to have a chance at existence in the universe, we need to figure out either how to keep this superorganism alive or help it reproduce.

We each have to recognize our own individual limitations of time, commitment, and brainpower. So, I’m trying to limit my work to the study of biological evolution rather than conjuring up a ‘pet theory of everything’. However, I’m also trying not to let those disciplinary and institutional barriers limit the tools I find valuable, or the people I interact with.

So, the more I’ve thought about the complexity of evolution (for now let’s just say ‘complexity’ = ‘anything humans don’t yet understand’), the more I’ve been driven to search for new languages. And in that search, I’ve been driven toward pure mathematics, where there are many exciting languages lurking around. Perhaps one of these languages has already obviated the need to invent new ideas to understand biological evolution… or perhaps an altogether new language needs to be constructed.

The prospects of a general theory of evolution point to the same intellectual challenge that we see in the quote above from Bar-Yam: assuming we’d like to be able to consistently manipulate the universe, when can we neglect details and when can’t we?

Consider the level of organization concept. Since different details of a system can be effectively ignored at different scales, our scientific theories have themselves become ‘stratified’:

• G. L. Farre, The energetic structure of observation: a philosophical disquisition, American Behavioral Scientist 40 (May 1997), 717-728.

In other words, science tends to be organized in ‘layers’. These layers have come to be conceived of as levels of organization, and each scientific theory tends to address only one of these levels (click the image to see the flash animation that ascends through many scales or levels):

It might be useful to work explicitly on connecting theories that tell us about particular levels of organization in order to attempt to develop some theories that transcend levels of organization. One type of insight that could be gained from this approach is an understanding of the mutual development of bottom-up ostensibly mechanistic models of simple systems and top-down initially phenomenological models of complex ones.

Simon has written an interesting discussion of the quasi-continuum that ranges from simple systems to complex ones:

• H. A. Simon, The architecture of complexity, Proceedings of the American Philosophical Society 106 (1962), 467–482.

But if we take an ideological perspective on science that says “let’s unify everything!” (scientific monism), a significant challenge is the development of a language able to unify our descriptions of simple and complex systems. Such a language might help communication among scientists who work with complex systems that apparently involve multiple levels of organization. Something like category theory may provide the nucleus of the framework necessary to formally address this challenge. But, in order to head in that direction, I’ll try out a few examples in a series of posts, albeit from the somewhat limited perspective of a biologist, from which some patterns might begin to surface.

In this introductory post, I’ll try to set a basis for thinking about this tension between simple and complex systems without wading through any treatises on ‘complexity’. It will be remarkably imprecise, but I’ll try to describe the ways in which I think it provides a useful metaphor for thinking about how we humans have dealt with this simple ↔ complex tension in science.

Another tack that I think could accomplish a similar goal, perhaps in a clearer way, would be to discuss fractals, power laws and maybe even renormalization. I might try that out in a later post if I get a little help from my new Azimuth friends, but I don’t think I’m qualified yet to do it alone.

Simple and complex systems

What is the organizational structure of the products of evolutionary processes? Herbert Simon provides a perspective that I find intuitive in his parable of two watchmakers.

He argues that the systems containing modules that don’t instantaneously fall apart (‘stable intermediates’) and can be assembled hierarchically take less time to evolve complexity than systems that lack stable intermediates. Given a particular set of internal and environmental constraints that can only be satisfied by some relatively complex system, a hierarchically organized one will be capable of meeting those constraints with the fewest resources and in the least time (i.e. most efficiently). The constraints any system is subject to determine the types of structures that can form. If hierarchical organization is an unavoidable outcome of evolutionary processes, it should be possible to characterize the causes that lead to its emergence.

Simon describes a property that some complex systems have in common, which he refers to as ‘near decomposability’:

• H. A. Simon, Near decomposability and the speed of evolution, Industrial and Corporate Change 11 (June 2002), 587-599.

A system is nearly decomposable if it’s made of parts that interact rather weakly with each other; these parts in turn being made of smaller parts with the same property.

For example, suppose we have a system modelled by a first-order linear differential equation. To be concrete, consider the fictitious building imagined by Simon: the Mellon Institute, with 12 rooms. Suppose the temperature of the ith room at time t is T_i(t). Of course most real systems seem to be nonlinear, but for the sake of this metaphor we can imagine that the temperatures of these rooms interact in a linear way, like this:

\displaystyle{ \frac{d}{d t}T_i(t) = \sum_{j}a_{ij}\left(T_{j}(t)-T_{i}(t)\right),}

where a_{ij} are some numbers. Suppose also that the matrix a_{ij} looks like this:

For the sake of the metaphor I’m trudging through here, let’s also assume

a\;\gg\;\epsilon_l\;\gg\;\epsilon_2

Then our system is nearly decomposable. Why? It has three ‘layers’, with two cells at the top level, each divided into two subcells, and each of these subdivided into three sub-subcells. The numbers of the rows and columns designate the cells, cells 1–6 and 7–12 constitute the two top-level subsystems, cells 1–3, 4–6, 7–9 and 10–12 the four second-level sub- systems. The interactions within the latter subsystems have intensity a, those within the former two subsystems, intensity \epsilon_l, and those between components of the largest subsystems, intensity \epsilon_2. This is why Simon states that this matrix is in near-diagonal form. Another, probably more common, terminology for this would be near block diagonal form. This terminology is a bit sloppy, but it basically means that we have a square matrix whose diagonal entries are square matrices and all other elements are approximately zero. That ‘approximately’ there is what differentiates near block diagonal matrices from honest block diagonal matrices whose off diagonal matrix elements are precisely zero.

This is a trivial system, but it illustrates that the near decomposability of the coefficient matrix allows these equations to be solved in a near hierarchical fashion. As an approximation, rather than simulating all the equations at once (e.g. all twelve in this example) one can take a recursive approach and solve the four systems of three equations (each of the blocks containing a‘s), and average the results to produce initial conditions for two systems of two equations with coefficients:

\begin{array}{cccc} \epsilon_1 & \epsilon_1 & \epsilon_2 & \epsilon_2 \\         \epsilon_1 & \epsilon_1 & \epsilon_2 & \epsilon_2\\     \epsilon_2 & \epsilon_2 & \epsilon_1 & \epsilon_1 \\     \epsilon_2 & \epsilon_2 & \epsilon_1 & \epsilon_1        \end{array}

and then average those results to produce initial conditions for a single system of two equations with coefficients:

\begin{array}{cc} \epsilon_2 & \epsilon_2 \\     \epsilon_2 & \epsilon_2    \end{array}

This example of simplification indicates that the study of a nearly decomposable systems system can be reduced to a series of smaller modules, which can be simulated in less computational time, if the error introduced in this approximation is tolerable. The degree to which this method saves time depends on the relationship between the size of the whole system and the size and number of hierarchical levels. However, as an example, given that the time complexity for matrix inversion (i.e. solving a system of linear equations) is O(n^2), then the hierarchical decomposition would lead to an algorithm with time complexity

\displaystyle{ O\left(\left(\frac{n}{L}\right)^2\right)}

where L is the number of levels in the decomposition. (For example, L=4 in the Mellon Institute, assuming the individual rooms are the lowest level).

All of this deserves to be made much more precise. However, there are some potential metaphorical consequences for the evolution of complex systems:

If we begin with a population of systems of comparable complexity, some of which are nearly decomposable and some of which are not, the nearly decomposable systems will, on average, increase their fitness through evolutionary processes much faster than the remaining systems, and will soon come to dominate the entire population. Notice that the claim is not that more complex systems will evolve more rapidly than less complex systems, but that, at any level of complexity, nearly decomposable systems will evolve much faster than systems of comparable complexity that are not nearly decomposable. – Herbert Simon, 2002

The point I’d like to make is that in this system, the idea of switching back and forth between simple and complex perspectives is made explicit: we get a sort of conceptual parallax:

In this simple case, the approximation that Simon suggests works well; however, for some other systems, it wouldn’t work at all. If we aren’t careful, we might even become victims of the Dunning-Kruger effect. In other words: if we don’t understand a system well from the start, we may overestimate how well we understand the limitations inherent to the simplifications we employ in studying it.

But if we at least recognize the potential of falling victim to the Dunning-Kruger effect, we can vigilantly guard against it in trying to understand, for example, the currently paradoxical tension between ‘groups’ and ‘individuals’ that lies at the heart of evolutionary theory… and probably also the caricatures of evolution that breed social controversy.

Keeping this in mind, my starting point in the next post in this series will be to provide some examples of hierarchical organization in biological systems. I’ll also set the stage for a discussion of evolution viewed as a dynamic process involving structural and functional transitions in hierarchical organization—or for the physicists out there, something like phase transitions!


Eddy Who?

25 August, 2011

Or: A very short introduction to turbulence

guest post by Tim van Beek

Have a look at this picture:

Then look at this one:

Do they look similar?

They should! They are both examples of a Kelvin-Helmoltz instability.

The first graphic is a picture of billow clouds (the fancier name is altostratus undulatus clouds):

The picture is taken from:

• C. Donald Ahrens: Meteorology Today, 9th edition, Brooks/Cole, Kentucky, 2009.

The second graphic:

shows a lab experiment and is taken from:

• G.L. Brown and A. Roshko, online available Density effects and large structure in turbulent mixing layers, Journal of Fluid Mechanics 64 (1974), 775-816.

Isn’t it strange that clouds in the sky would show the same pattern as some gases in a small laboratory experiment? The reason for this is not quite understood today. In this post, I would like to talk a little bit about what is known.

Matter that tries to get out of its own way

Fluids like water and air can be well described by Newton’s laws of classical mechanics. When you start learning about classical mechanics, you consider discrete masses, most of the time. Billiard balls, for example. But it is possible to formulate Newton’s laws of motion for fluids by treating them as ‘infinitely many infinitesimally small’ billiard balls, all pushing and rubbing against each other and therefore trying to get out of the way of each other.

If we do this, we get some equations describing fluid flow: the Navier-Stokes equations.

The Navier-Stokes equations are a complicated set of nonlinear partial differential equations. A lot of mathematical questions about them are still unanswered, like: under what conditions is there a smooth solution to these equations? If you can answer that question, you will win one of the a million dollars from the Clay Mathematics Institute.

If you completed the standard curriculum of physics as I did, chances are that you never attended a class on fluid dynamics. At least I never did. When you take a first look at the field, you will notice: the literature about the Navier-Stokes equations alone is huge! Not to mention all the special aspects of numerical simulations, special aspects of the climate and so on.

So it is nice to find a pedagogical introduction to the subject for people who have some background knowledge in partial differential equations, for the mathematical theory:

• C. Foias, R. Rosa, O. Manley and R. Temam, Navier-Stokes Equations and Turbulence, Cambridge U. Press, Cambridge, 2001.

So, there is a lot of fun to be had for the mathematically inclined. But today I would like to talk about an aspect of fluid flows that also has a tremendous practical importance, especially for the climate of the Earth: turbulence!

There is no precise definition of turbulence, but people know it when they see it. A fluid can flow in layers, with one layer above the other, maybe slightly slower or faster. Material of one layer does hardly mix with material of another layer. These flows are called laminar flows. When a laminar flow gets faster and faster, it turns into a turbulent flow at some point:

This is a fluid flow inside a circular pipe, with a layer of some darker fluid in the middle.

As a first guess we could say that a characteristic property of turbulent flow is the presence of circular flows, commonly called eddies.

Tempests in Teapots

A funny aspect of the Navier-Stokes equations is that they don’t come with any recommended length scale. Properties of the fluid flow like velocity and pressure are modelled as smooth functions of continuous time and space. Of course we know that this model does not work on a atomic length scale, where we have to consider individual atoms. Pressure and velocity of a fluid flow don’t make any sense on a length scale that is smaller than the average distance between electrons and the atom nucleus.

We know this, but this is a fact that is not present in the model comprised by the Navier-Stokes equations!

But let us look at bigger length scales. An interesting feature of the solutions of the Navier-Stokes equations is that there are fluid flows that stretch over hundreds of meters that look like fluid flows that stretch over centimeters only. And it is really astonishing that this phenomenon can be observed in nature. This is another example of the unreasonable effectiveness of mathematical models.

You have seen an example of this in the introduction already. That was a boundary layer instability. Here is a full blown turbulent example:

The last two pictures are from the book:

• Arkady Tsinober: An Informal Conceptual Introduction to Turbulence, 2nd edition, Springer, Fluid Mechanics and Its Applications Volume 92, Berlin, 2009.

This is a nice introduction to the subject, especially if you are more interested in phenomenology than mathematical details.

Maybe you noticed the “Reynolds number” in the label text of the last picture. What is that?

People in business administration like management ratios; they throw all the confusing information they have about a company into a big mixer and extract one or two numbers that tell them where they stand, like business volume and earnings. People in hydrodynamics are somewhat similar; they define all kinds of “numbers” that condense a lot of information about fluid flows.

A CEO would want to know if the earnings of his company are positive or negative. We would like to know a number that tells us if a fluid flow is laminar or turbulent. Luckily, such a number already exists. It is the Reynolds number! A low number indicates a laminar flow, a high number a turbulent flow. Like the calculation of the revenue of a company, the calculation of the Reynolds number of a given fluid flow is not an exact science. Instead there is some measure of estimation necessary. The definition involves, for example, a “characteristic length scale”. This is a fuzzy concept that usually involves some object that interacts with – in our case – the fluid flow. The characteristic length scale in this case is the physical dimension of the object. While there is usually no objectively correct way to assing a “characteristic length” to a three dimensional object, this concept allows us nevertheless to distinguish the scattering of water waves on a ocean liner (length scale ≈ 103 meter) from their scattering on a peanut (length scale ≈ 10-2 meter).

The following graphic shows laminar and turbulent flows and their characteristic Reynolds numbers:

schematic display of laminar and turbulent flow for different Reynolds numbers

This graphic is from the book

• Thomas Bohr, Mogens H. Jensen, Giovanni Paladin and Angelo Vulpiani, Dynamical Systems Approach to Turbulence, Cambridge U. Press, Cambridge, 1998

But let us leave the Reynolds number for now and turn to one of its ingredients: viscosity. Understanding viscosity is important for understanding how eddies in a fluid flow are connected to energy dissipation.

Eddies dissipate kinetic energy

“Eddies,” said Ford, “in the space-time continuum.”

“Ah,” nodded Arthur, “is he? Is he?” He pushed his hands into the pocket of his dressing gown and looked knowledgeably into the distance.

“What?” said Ford.

“Er, who,” said Arthur, “is Eddy, then, exactly, then?”

– from Douglas Adams: Life, the Universe and Everything

A fluid flow can be pictured as consisting of a lot of small fluid packages that move alongside each other. In many situations, there will be some friction between these packages. In the case of fluids, this friction is called viscosity.

It is an empirical fact that at small velocities fluid flows are laminar: there are layers upon layers, with one layer moving at a constant speed, and almost no mixing. At the boundaries, the fluid will attach to the surrounding material, and the relative fluid flow will be zero. If you picture such a flow between a plate that is at rest, and a plate that is moving forward, you will see that due to friction between the layers a force needs to be exerted to keep the moving plate moving:

In the simplest approximation, you will have to exert some force F per unit area A, in order to sustain a linear increase of the velocity of the upper plate along the y-axis, \partial u / \partial y. The constant of proportionality is called the viscosity \mu:

\displaystyle{ \frac{F}{A} = \mu \frac{\partial u}{\partial y}}

More friction means a bigger viscosity: honey has a bigger viscosity than water.

If you stir honey, the fluid flow will come to a halt rather fast. The energy that you put in to start the fluid flow is turned to heat by dissipation. This mechanism is of course related to friction and therefore to viscosity.

It is possible to formulate an exact formula for this dissipation process using the Navier-Stokes equations. It is not hard to prove it, but I will only explain the involved gadgets.

A fluid flow in three dimensions can be described by stating the velocity of the fluid flow at a certain time t and \vec{x} \in \mathbb{R}^3 (I’m not specifying the region of the fluid flow or any boundary or initial conditions). Let’s call the velocity \vec{u}(t, \vec{x}).

Let’s assume that the fluid has a constant density \rho. Such a fluid is called incompressible. For convenience let us assume that the density is 1: \rho = 1. Then the kinetic energy E(\vec{u}) of a fluid flow at a fixed time t is given by

\displaystyle{ E(\vec{u}) = \int \| \vec{u}(t, \vec{x}) \|^2 \; d^3 x }

Let’s just assume that this integral is finite for the moment. This is the first gadget we need.

The second gadget is called enstrophy \epsilon of the fluid flow. This is a measure of how much eddies there are. It is the integral

\displaystyle{\epsilon = \int \| \nabla \times \vec{u} \|^2 \; d^3 x }

where \nabla \times denotes the curl of the fluid velocity. The faster the fluid rotates, the bigger the curl is.

(The math geeks will notice that the vector fields \vec{u} that have a finite kinetic energy and a finite enstrophy are precisely the elements of the Sobolev space H^1(\mathbb{R}^3))

Here is the relationship of the decay of the kinetic energy and the enstrophy, which is a consequence of the Navier-Stokes equations (and suitable boundary conditions):

\displaystyle{\frac{d}{d t} E = - \mu \epsilon}

This equation says that the energy decays with time, and it decays faster if there is a higher viscosity, and if there are more and stronger eddies.

If you are interested in the mathematically precise derivation of this equation, you can look it up in the book I already mentioned:

• C. Foias, R. Rosa, O. Manley and R. Temam, Navier-Stokes Equations and Turbulence, Cambridge U. Press, Cambridge, 2001.

This connection of eddies and dissipation could indicate that there is also a connection of eddies and some maximum entropy principle. Since eddies maximize dissipation, natural fluid flows should somehow tend towards the production of eddies. It would be interesting to know more about this!

In this post we have seen eddies at different length scales. There are buzzwords in meteorology for this:

You have seen eddies at the ‘microscale’ (left) and at the ‘mesoscale’ (middle). A blog post about eddies should of course mention the most famous eddy of the last decade, which formed at the ‘synoptic scale’:

Do you recognize it? That was Hurricane Katrina.

It is obviously important to understand disasters like this one on the synoptic scale. This is an active topic of ongoing research, both in meteorology and in climate science.


This Week’s Finds (Week 318)

21 August, 2011

The Earth’s climate is complicated. How well do we understand how it will react to the drastic increase in atmospheric carbon dioxide that we’re imposing on it? One reason it’s hard to be 100% sure is that the truly dramatic changes most scientists expect lie mostly in the future. There’s a lot of important evidence we’ll get only when it’s too late.

Luckily, the Earth’s past also shows signs of dramatic climate change: for example, the glacial cycles I began discussing last time in "week317". These cycles make an interesting test of how well we understand climate change. Of course, their mechanism is very different from that of human-caused global warming, so we might understand one but not the other. Indeed, earlier episodes like the Paleocene-Eocene Thermal Maximum might shed more light on what we’re doing to the Earth now! But still, the glacial cycles are an impressive instance of dramatic climate change, which we’d do well to understand.

As I hinted last week, a lot of scientists believe that the Earth’s glacial cycles are related to cyclic changes in the Earth’s orbit and the tilt of its axis. Since one of the first scientists to carefully study this issue was Milutin Milankovitch, these are called Milankovitch cycles. The three major types of Milankovitch cycle are:

• changes in the eccentricity of the Earth’s orbit – that is, how much the orbit deviates from being a circle:




(changes greatly exaggerated)

• changes in the obliquity, or tilt of the Earth’s axis:



precession, meaning changes in the direction of the Earth’s axis relative to the fixed stars:



Now, the first important thing to realize is this: it’s not obvious that Milankovitch cycles can cause glacial cycles. During a glacial period, the Earth is about 5°C cooler than it is now. But the Milankovitch cycles barely affect the overall annual amount of solar radiation hitting the Earth!

This fact is clear for precession or changes in obliquity, since these just involve the tilt of the Earth’s axis, and the Earth is nearly a sphere. The amount of Sun hitting a sphere doesn’t depend on how the sphere is ’tilted’.

For changes in the eccentricity of the Earth’s orbit, this fact is a bit less obvious. After all, when the orbit is more eccentric, the Earth gets closer to the Sun sometimes, but farther at other times. So you need to actually sit down and do some math to figure out the net effect. Luckily, Greg Egan did this for us—I’ll show you his calculation at the end of this article. It turns out that when the Earth’s orbit is at its most eccentric, it gets very, very slightly more energy from the Sun each year: 0.167% more than when its orbit is at its least eccentric.

So, there are interesting puzzles involved in the Milankovitch cycles. They don’t affect the total amount of radiation that hits the Earth each year—not much, anyway—but they do cause substantial changes in the amount of radiation that hits the Earth at various different latitudes in various different seasons. We need to understand what such changes might do.

James Croll was one of the first to think about this, back around 1875. He decided that what really matters is the amount of sunlight hitting the far northern latitudes in winter. When this was low, he claimed, glaciers would tend to form and an ice age would start. But later, in the 1920s, Milankovitch made the opposite claim: what really matters is the amount of sunlight hitting the far northern latitudes in summer. When this was low, an ice age would start.

If we take a quick look at the data, we see that the truth is not obvious:


I like this graph because it’s pretty… but I wish the vertical axes were labelled. We will see some more precise graphs in future weeks.

Nonetheless, this graph gives some idea of what’s going on. Precession, obliquity and eccentricity vary in complex but still predictable ways. From this you can compute the amount of solar energy that hits the surface of the Earth’s atmosphere on July 1st at a latitude of 65° N. That’s the yellow curve. People believe this quantity has some relation to the Earth’s temperature, as shown by the black curve at bottom. However, the relation is far from clear!

Indeed, if you only look at this graph, you might easily decide that Milankovitch cycles are not important in causing glacial cycles. But people have analyzed temperature proxies over long spans of time, and found evidence for cyclic changes at periods that match those of the Milankovitch cycles. Here’s a classic paper on this subject:

• J. D. Hays, J. Imbrie, and N. J. Shackleton, Variations in the earth’s orbit: pacemaker of the Ice Ages, Science 194 (1976), 1121-1132.

They selected two sediment cores from the Indian ocean, which contain sediments deposited over the last 450,000 years. They measured:

1) Ts, an estimate of summer sea-surface temperatures at the core site, derived from a statistical analysis of tiny organisms called radiolarians found in the sediments.

2) δ18O, the excess of the heavy isotope of oxygen in tiny organisms called foraminifera also found in the sediments.

3) The percentage of radiolarians that are Cycladophora davisiana—a certain species not used in the estimation of Ts.

Identical samples were analyzed for the three variables at 10-centimeter intervals throughout each core. Then they took a Fourier transform of this data to see at which frequencies these variables wiggle the most! When we take the Fourier transform of a function and then square it, the result is called the power spectrum. So, they actually graphed the power spectra for these three variables:

The top graph shows the power spectra for Ts, δ18O, and the percentage of Cycladophora davisiana. The second one shows the spectra after a bit of extra messing around. Either way, there seem to be peaks at frequencies of 19, 23, 42 and roughly 100 thousand years. However the last number is quite fuzzy: if you look, you’ll see the three different power spectra have peaks at 94, 106 and 122 thousand years.

So, some sort of cycles seem to be occurring. This is far from the only piece of evidence, but it’s a famous one.

Now let’s go over the three major forms of Milankovitch cycle, and keep our eye out for cycles that take place every 19, 23, 42 or roughly 100 thousand years!

Eccentricity

The Earth’s orbit is an ellipse, and the eccentricity of this ellipse says how far it is from being circular. But the eccentricity of the Earth’s orbit slowly changes: it varies from being nearly circular, with an eccentricity of 0.005, to being more strongly elliptical, with an eccentricity of 0.058. The mean eccentricity is 0.028. There are several periodic components to these variations. The strongest occurs with a period of 413,000 years, and changes the eccentricity by ±0.012. Two other components have periods of 95,000 and 123,000 years.

The eccentricity affects the percentage difference in incoming solar radiation between the perihelion, the point where the Earth is closest to the Sun, and the aphelion, when it is farthest from the Sun. This works as follows. The percentage difference between the Earth’s distance from the Sun at perihelion and aphelion is twice the eccentricity, and the percentage change in incoming solar radiation is about twice that. The first fact follows from the definition of eccentricity, while the second follows from differentiating the inverse-square relationship between brightness and distance.

Right now the eccentricity is 0.0167, or 1.67%. Thus, the distance from the Earth to Sun varies 3.34% over the course of a year. This in turn gives an annual variation in incoming solar radiation of about 6.68%. Note that this is not the cause of the seasons: those arise due to the Earth’s tilt, and occur at different times in the northern and southern hemispheres.

Obliquity

The angle of the Earth’s axial tilt with respect to the plane of its orbit, called the obliquity, varies between 22.1° and 24.5° in a roughly periodic way, with a period of 41,000 years. When the obliquity is high, the strength of seasonal variations is stronger.

Right now the obliquity is 23.44°, roughly halfway between its extreme values. It is decreasing, and will reach its minimum value around the year 10,000 CE.

Precession

The slow turning in the direction of the Earth’s axis of rotation relative to the fixed stars, called precession, has a period of roughly 23,000 years. As precession occurs, the seasons drift in and out of phase with the perihelion and aphelion of the Earth’s orbit.

Right now the perihelion occurs during the southern hemisphere’s summer, while the aphelion is reached during the southern winter. This tends to make the southern hemisphere seasons more extreme than the northern hemisphere seasons.

The gradual precession of the Earth is not due to the same physical mechanism as the wobbling of the top. That sort of wobbling does occur, but it has a period of only 427 days. The 23,000-year precession is due to tidal interactions between the Earth, Sun and Moon. For details, see:

• John Baez, The wobbling of the Earth and other curiosities.

In the real world, most things get more complicated the more carefully you look at them. For example, precession actually has several periodic components. According to André Berger, a top expert on changes in the Earth’s orbit, the four biggest components have these periods:

• 23,700 years

• 22,400 years

• 18,980 years

• 19,160 years

in order of decreasing strength. But in geology, these tend to show up either as a single peak around the mean value of 21,000 years, or two peaks at frequencies of 23,000 and 19,000 years.

To add to the fun, the three effects I’ve listed—changes in eccentricity, changes in obliquity, and precession—are not independent. According to Berger, cycles in eccentricity arise from ‘beats’ between different precession cycles:

• The 95,000-year eccentricity cycle arises from a beat between the 23,700-year and 19,000-year precession cycles.

• The 123,000-year eccentricity cycle arises from a beat between the 22,4000-year and 18,000-year precession cycles.

We should delve into all this stuff more deeply someday. For now, let me just refer you to this classic review paper:

• André Berger, Pleistocene climatic variability at astronomical frequencies, Quaternary International 2 (1989), 1-14.

Later, as I get up to speed, I’ll talk about more modern work.

Paleontology versus astronomy

So now we can compare the data from ocean sediments to the Milankovitch cycles as computed in astronomy:

• The roughly 19,000-year cycle in ocean sediments may come from 18,980-year and 19,160-year precession cycles.

• The roughly 23,000-year cycle in ocean sediments may come from 23,700-year precession cycle.

• The roughly 42,000-year cycle in ocean sediments may come from the 41,000-year obliquity cycle.

• The roughly 100,000-year cycle in ocean sediments may come from the 95,000-year and 123,000-year eccentricity cycles.

Again, the last one looks the most fuzzy. As we saw, different kinds of sediments seem to indicate cycles of 94, 106 and 122 thousand years. At least two of these periods match eccentricity cycles fairly well. But a detailed analysis would be required to distinguish between real effects and coincidences in this subject!

The effect of eccentricity

I bet some of you are hungry for some actual math. As I mentioned, it takes some work to see how changes in the eccentricity of the Earth’s orbit affect the annual average of sunlight hitting the top of the Earth’s atmosphere. Luckily Greg Egan has done this work for us. While the result is surely not new, his approach makes nice use of the fact that both gravity and solar radiation obey an inverse-square law. That’s pretty cool.

Here is his calculation:

The angular velocity of a planet is

\displaystyle{\frac{d \theta}{d t} = \frac{J}{m r^2} }

where J is the constant orbital angular momentum of the planet and m is its mass. Thus the radiant energy delivered per unit time to the planet is

\displaystyle{ \frac{d U}{d t} = \frac{C}{r^2}}

for some constant C. It follows that the energy delivered per unit of angular progress around the orbit is

\displaystyle{ \frac{d U}{d \theta} = \frac{C}{r^2} \frac{d t}{d \theta} = \frac{C m}{J} }

So, the total energy delivered in one period will be

\displaystyle{ U=\frac{2\pi C m}{J} }

How can we relate the orbital angular momentum J to the shape of the orbit? If you equate the total energy of the planet, kinetic \frac{1}{2}m v^2 plus potential -\frac{G M m}{r}, at its aphelion r_1 and perihelion r_2, and use J to get the velocity in the kinetic energy term from its distance, v=\frac{J}{m r}, when we solve for J we get:

\displaystyle{J = m \sqrt{\frac{2 G M r_1 r_2}{r_1+r_2}} = m b \sqrt{\frac{G M}{a}}}

where

\displaystyle{ a=\frac{1}{2} (r_1+r_2)}

is the semi-major axis of the orbit and

\displaystyle{ b=\sqrt{r_1 r_2} }

is the semi-minor axis. But we can also relate J to the period of the orbit, T, by integrating the rate at which orbital area is swept out by the planet:

\displaystyle{\frac{1}{2}  r^2 \frac{d \theta}{d t} = \frac{J}{2 m} }

over one orbit. Since the area of an ellipse is \pi a b, this gives us:

\displaystyle{ J = \frac{2 \pi a b m}{T} }

Equating these two expressions for J shows that the period is:

\displaystyle{ T = 2 \pi \sqrt{\frac{a^3}{G M}}}

So the period depends only on the semi-major axis; for a fixed value of a, it’s independent of the eccentricity.

As the eccentricity of the Earth’s orbit changes, the orbital period T, and hence the semi-major axis a, remains almost constant. So, we have:

\displaystyle{ U=\frac{2\pi C m}{J} = \frac{2\pi C}{b} \sqrt{\frac{a}{G M}} }

Expressing the semi-minor axis in terms of the semi-major axis and the eccentricity, b^2 = a^2 (1-e^2), we get:

\displaystyle{U=\frac{2\pi C}{\sqrt{G M a (1-e^2)}}}

So to second order in e, we have:

\displaystyle{U = \frac{\pi C}{\sqrt{G M a}} (2+e^2) }

The expressions simplify if we consider average rate of energy delivery over an orbit, which makes all the grungy constants related to gravitational dynamics go away:

\displaystyle{\frac{U}{T} = \frac{C}{a^2 \sqrt{1-e^2}} }

or to second order in e:

\displaystyle{ \frac{U}{T} = \frac{C}{a^2} (1+\frac{1}{2} e^2) }

We can now work out how much the actual changes in the Earth’s orbit affect the amount of solar radiation it gets! The eccentricity of the Earth’s orbit varies between 0.005 and 0.058. The total energy the Earth gets each year from solar radiation is proportional to

\displaystyle{ \frac{1}{\sqrt{1-e^2}} }

where e is the eccentricity. When the eccentricity is at its lowest value, e = 0.005, we get

\displaystyle{ \frac{1}{\sqrt{1-e^2}} = 1.0000125 }

When the eccentricity is at its highest value, e = 0.058, we get

\displaystyle{\frac{1}{\sqrt{1-e^2}} = 1.00168626 }

So, the change is about

\displaystyle{1.00168626/1.0000125 = 1.00167373 }

In other words, a change of merely 0.167%.

That’s very small And the effect on the Earth’s temperature would naively be even less!

Naively, we can treat the Earth as a greybody: an ideal object whose tendency to absorb or emit radiation is the same at all wavelengths and temperatures. Since the temperature of a greybody is proportional to the fourth root of the power it receives, a 0.167% change in solar energy received per year corresponds to a percentage change in temperature roughly one fourth as big. That’s a 0.042% change in temperature. If we imagine starting with an Earth like ours, with an average temperature of roughly 290 kelvin, that’s a change of just 0.12 kelvin!

The upshot seems to be this: in a naive model without any amplifying effects, changes in the eccentricity of the Earth’s orbit would cause temperature changes of just 0.12 °C!

This is much less than the roughly 5 °C change we see between glacial and interglacial periods. So, if changes in eccentricity are important in glacial cycles, we have some explaining to do. Possible explanations include season-dependent phenomena and climate feedback effects. Probably both are very important!

Next time I’ll start talking about some theories of how Milankovitch cycles might cause the glacial cycles. I thank Frederik De Roo, Martin Gisser and Cameron Smith for suggesting improvements to this issue before its release, over on the Azimuth Forum. Please join us over there.


Little did I suspect, at the time I made this resolution, that it would become a path so entangled that fully twenty years would elapse before I could get out of it. – James Croll, on his decision to study the cause of the glacial cycles


Bayesian Computations of Expected Utility

19 August, 2011

GiveWell is an organization that rates charities. They’ve met people who argue that

charities working on reducing the risk of sudden human extinction must be the best ones to support, since the value of saving the human race is so high that “any imaginable probability of success” would lead to a higher expected value for these charities than for others.

For example, say I have a dollar to spend on charity. One charity says that with this dollar they can save the life of one child in Somalia. Another says that with this dollar they can increase by .000001% our chance of saving 1 billion people from the effects of a massive asteroid colliding with the Earth.

Naively, in terms of the expected number of lives saved, the latter course of action seems 10 times better, since

.000001% × 1 billion = 10

But is it really better?

It’s a subtle question, with all sorts of complicating factors, like why should I trust these guys?

I’m not ready to present a thorough analysis of this sort of question today. But I would like to hear what you think about it. And I’d like you to read what the founder of Givewell has to say about it:

• Holden Karnofsky, Why we can’t take expected value estimates literally (even when they’re unbiased), 18 August 2011.

He argues against what he calls an ‘explicit expected value’ or ‘EEV’ approach:

The mistake (we believe) is estimating the “expected value” of a donation (or other action) based solely on a fully explicit, quantified formula, many of whose inputs are guesses or very rough estimates. We believe that any estimate along these lines needs to be adjusted using a “Bayesian prior”; that this adjustment can rarely be made (reasonably) using an explicit, formal calculation; and that most attempts to do the latter, even when they seem to be making very conservative downward adjustments to the expected value of an opportunity, are not making nearly large enough downward adjustments to be consistent with the proper Bayesian approach.

His focus, in short, is on the fact that anyone saying “this money can increase by .000001% our chance of saving 1 billion people from an asteroid impact” is likely to be pulling those numbers from thin air. If they can’t really back up their numbers with a lot of hard evidence, then our lack of confidence in their estimate should be taken into account somehow.

His article spends a lot of time analyzing a less complex but still very interesting example:

It seems fairly clear that a restaurant with 200 Yelp reviews, averaging 4.75 stars, ought to outrank a restaurant with 3 Yelp reviews, averaging 5 stars. Yet this ranking can’t be justified in an explicit expected utility framework, in which options are ranked by their estimated average/expected value.

This is the only question I really want to talk about today. Actually I’ll focus on a similar question that Tim van Beek posed on this blog:

You have two kinds of fertilizer, A and B. You know that of 4 trees who got A, three thrived and one died. Of 36 trees that got B, 24 thrived and 12 died. Which fertilizer would you buy?

So, 3/4 of the trees getting fertilizer A thrived, while only 2/3 of those getting fertilizer B thrived. That makes fertilizer A seem better. However, the sample size is considerably larger for fertilizer B, so we may feel more confident about the results in this case. Which should we choose?

Nathan Urban tackled the problem in an interesting way. Let me sketch what he did before showing you his detailed work.

Suppose that before doing any experiments at all, we assume the probability \pi that a fertilizer will make a tree thrive is a number uniformly distributed between 0 and 1. This assumption is our “Bayesian prior”.

Note: I’m not saying this prior is “correct”. You are allowed to choose a different prior! Choosing a different prior will change your answer to this puzzle. That can’t be helped. We need to make some assumption to answer this kind of puzzle; we are simply making it explicit here.

Starting from this prior, Nathan works out the probability that \pi has some value given that when we apply the fertilizer to 4 trees, 3 thrive. That’s the black curve below. He also works out the probability that \pi has some value given that when we apply the fertilizer to 36 trees, 24 thrive. That’s the red curve:

The red curve, corresponding to the experiment with 36 trees, is much more sharply peaked. That makes sense. It means that when we do more experiments, we become more confident that we know what’s going on.

We still have to choose a criterion to decide which fertilizer is best! This is where ‘decision theory’ comes in. For example, suppose we want to maximize the expected number of the trees that thrive. Then Nathan shows that fertilizer A is slightly better, despite the smaller sample size.

However, he also shows that if fertilizer A succeeded 4 out of 5 times, while fertilizer B succeeded 7 out of 9 times, the same evaluation procedure would declare fertilizer B better! Its percentage success rate is less: about 78% instead of 80%. However, the sample size is larger. And in this particular case, given our particular Bayesian prior and given what we are trying to maximize, that’s enough to make fertilizer B win.

So if someone is trying to get you to contribute to a charity, there are many interesting issues involved in deciding if their arguments are valid or just a bunch of… fertilizer.

Here is Nathan’s detailed calculation:

It’s fun to work out an official ‘correct’ answer mathematically, as John suggested. Of course, this ends up being a long way of confirming the obvious—and the answer is only as good as the assumptions—but I think it’s interesting anyway. In this case, I’ll work it out by maximizing expected utility in Bayesian decision theory, for one choice of utility function. This dodges the whole risk aversion point, but it opens discussion for how the assumptions might be modified to account for more real-world considerations. Hopefully others can spot whether I’ve made mistakes in the derivations.

In Bayesian decision theory, the first thing you do is write down the data-generating process and then compute a posterior distribution for what is unknown.

In this case, we may assume the data-generating process (likelihood function) is a binomial distribution \mathrm{Bin}(s,n|\pi) for s successes in n trials, given a probability of success \pi. Fertilizer A corresponds to s=3, n=4 and fertilizer B corresponds to s=24, n=36.

The probability of success \pi is unknown, and we want to infer its posterior conditional on the data, p(\pi|s,n). To compute a posterior we need to assume a prior on \pi.

It turns out that the Beta distribution is conjugate to a binomial likelihood, meaning that if we assume a Beta distributed prior, the then the posterior is also Beta distributed. If the prior is \pi \sim \mathrm{Beta}(\alpha_0,\beta_0) then the posterior is

\pi \sim \mathrm{Beta}(\alpha=\alpha_0+s,\beta=\beta_0+n-s).

One choice for a prior is a uniform prior on [0,1], which corresponds to a \mathrm{Beta}(1,1) distribution. There are of course other prior choices which will lead to different conclusions. For this prior, the posterior is \mathrm{Beta}(\pi; s+1, n-s+1). The posterior mode is

(\alpha-1)/(\alpha+\beta-2) = s/n

and the posterior mean is

\alpha/(\alpha+\beta) = (s+1)/(n+2).

So, what is the inference for fertilizers A and B? I made a graph of the posterior distributions. You can see that the inference for fertilizer B is sharper, as expected, since there is more data. But the inference for fertilizer A tends towards higher success rates, which can be quantified.

Fertilizer A has a posterior mode of 3/4 = 0.75 and B has a mode of 2/3 = 0.667, corresponding to the sample proportions. The mode isn’t the only measure of central tendency we could use. The means are 0.667 for A and 0.658 for B; the medians are 0.686 for A and 0.661 for B. No matter which of the three statistics we choose, fertilizer A looks better than fertilizer B.

But we haven’t really done “decision theory” yet. We’ve just compared point estimators. Actually, we have done a little decision theory, implicitly. It turns out that picking the mean corresponds to the estimator which minimizes the expected squared error in \pi, where “squared error” can be thought of formally as a loss function in decision theory. Picking the median corresponds to minimizing the expected absolute loss, and picking the mode corresponds to minimizing the minimizing the 0-1 loss (where you lose nothing if you guess \pi exactly and lose 1 otherwise).

Still, these don’t really correspond to a decision theoretic view of the problem. We don’t care about the quantity \pi at all, let alone some point estimator of it. We only care about \pi indirectly, insofar as it helps us predict something about what the fertilizer will do to new trees. For that, we have to move from the posterior distribution p(\pi|s,n) to the predictive distribution

p(y|s,n) = \int p(y|\pi,n)\,p(\pi|s,n)\,d\pi ,

where y is a random variable indicating whether a new tree will thrive under treatment. Here I assume that the success of new trees follows the same binomial distribution as in the experimental group.

For a Beta posterior, the predictive distribution is beta-binomial, and the expected number of successes for a new tree is equal to the mean of the Beta distribution for \pi – i.e. the posterior mean we computed before, (s+1)/(n+2). If we introduce a utility function such that we are rewarded 1 util for a thriving tree and 0 utils for non-thriving tree, then the expected utility is equal to the expected success rate. Therefore, under these assumptions, we should choose the fertilizer that maximizes the quantity (s+1)/(n+2), which, as we’ve seen, favors fertilizer A (0.667) over fertilizer B (0.658).

An interesting mathematical question is, does this ever work out to a “non-obvious” conclusion? That is, if fertilizer A has a sample success rate which is greater than fertilizer B’s sample success rate, but expected utility maximization prefers fertilizer B? Mathematically, we’re looking for a set {s,s',n,n'} such that s/n>s'/n' but (s+1)/(n+2) < (s'+1)/(n'+2). (Also there are obvious constraints on s and s'.) The answer is yes. For example, if fertilizer A has 4 of 5 successes while fertilizer B has 7 of 9 successes.

By the way, on a quite different note: NASA currently rates the chances of the asteroid Apophis colliding with the Earth in 2036 at 4.3 × 10-6. It estimates that the energy of such a collision would be comparable with a 510-megatonne thermonuclear bomb. This is ten times larger than the largest bomb actually exploded, the Tsar Bomba. The Tsar Bomba, in turn, was ten times larger than all the explosives used in World War II.



There’s an interesting Chinese plan to deflect Apophis if that should prove necessary. It is, however, quite a sketchy plan. I expect people will make more detailed plans shortly before Apophis comes close to the Earth in 2029.


Follow

Get every new post delivered to your Inbox.

Join 3,095 other followers