A Toy Model of Climate Change

Introduction

A screenshot of the online book
Math for the People.

In Introductory Physics for Medicine and Biology, Russ Hobbie and I make use of toy models. Such mathematical descriptions are not intended to be accurate or realistic. Rather, they‘re simple models that capture the main idea without getting bogged down in the details. Today, I present an example of a toy model. It’s not related to medicine or biology, but instead describes climate change. I didn’t originally derive this model. Much of the analysis below comes from other sources, such as the online book Math for the People published by Mark Branson and Whitney George.

Earth Without an Atmosphere

First, consider the earth with no atmosphere. We will balance the energy coming into the earth from the sun with the energy from the earth that is radiated out into space. Our goal will be to calculate the earth’s temperature, T.

The power density (energy per unit time per unit area, in watts per square meter) emitted by the sun is called the solar constant, S. It depends on how far you are from the sun, but at the earth’s orbit S = 1360 W/m². To get the total power impinging on our planet, we must multiply S by the area subtended by the earth, which is πR², where R is the earth’s radius (R = 6.4 × 10⁶ meters). This gives SπR² = 1.8 × 10¹⁷ W, or nearly 200,000 TW (T, or tera-, means one trillion). That’s a lot of power. The total average power consumption by humanity is only about 20 TW, so there’s plenty of energy from the sun.

We often prefer to talk about the energy loss or gain per unit area of the earth’s surface. The surface area of the earth is 4πR² (the factor of four comes from the total surface area of the spherical earth, in contrast to the area subtended by the earth when viewed from the sun). The power per unit area of the earth’s surface is therefore SπR²/4πR², or S/4.

Not all of this energy is absorbed by the earth; some is reflected back into space. The albedo, a, is a dimensionless number that indicates the fraction of the sun’s energy that is reflected. The power absorbed per unit area is then (1 – a)S/4. About 30% of the sun’s energy is reflected (a = 0.3), so the power of sunlight absorbed by the earth per unit of surface area is 238 W/m².

What happens to that energy? The sun heats the earth to a temperature T. Any hot object radiates energy. Such thermal radiation is analyzed in Section 14.8 of Intermediate Physics for Medicine and Biology. The radiated power per unit area is equal to eσT⁴. The symbol σ is the Stefan-Boltzmann constant, σ = 5.7 × 10^–8 W/(m² K⁴). As stated earlier, T is the earth’s temperature. When raising the temperature to the fourth power, T must be expressed as the absolute temperature measured in kelvin (K). Sometimes it’s convenient at the end of a calculation to convert kelvin to the more familiar degrees Celsius (°C), where 0°C = 273 K. But remember, all calculations of T⁴ must use kelvin. Finally, e is the emissivity of the earth, which is a measure of how well the earth absorbs and emits radiation. The emissivity is another dimensionless number ranging between zero and one. The earth is an excellent emitter and absorber, so e = 1. From now on, I’ll not even bother including e in our equations, in which case the power density emitted is just σT⁴.

Let’s assume the earth is in steady state, meaning the temperature is not increasing or decreasing. Then the power in must equal the power out, so 

(1 – a)S/4 = σT⁴. 

Solving for the temperature gives

T = ∜[(1 – a)S/4σ] .

Because we know a, S, and σ, we can calculate the temperature. It is T = 254 K = –19°C. That’s really cold (remember, in the Celsius scale water freezes at 0°C). Without an atmosphere, the earth would be a frozen wasteland.

Earth With an Atmosphere

Often we can learn much from a toy model by adding in complications, one by one. Now, we’ll include an atmosphere around earth. We must keep track of the power into and out of both the earth and the atmosphere. The earth has temperature T_E and the atmosphere has temperature T_A.

First, let’s analyze the atmosphere. Sunlight passes right through the air without being absorbed because it’s mainly visible light and our atmosphere is transparent in the visible part of the spectrum. The main source of thermal (or infrared) radiation (for which the atmosphere is NOT transparent) is from the earth. We already know how much that is, σT_E⁴. The atmosphere only absorbs a fraction of the earth’s radiation, e_A, so the power per unit area absorbed by the atmosphere is e_AσT_E⁴.

Just like the earth, the atmosphere will heat up to a temperature T_A and emit its own thermal radiation. The emitted power per unit area is e_AσT_A⁴. However, the atmosphere has upper and lower surfaces, and we’ll assume they both emit equally well. So the total power emitted by the atmosphere per unit area is 2e_AσT_A⁴.

If we balance the power in and out of the atmosphere, we get 

e_AσT_E⁴ = 2e_AσT_A⁴. 

Interestingly, the fraction of radiation absorbed by the atmosphere, e_A, cancels out of our equation (a good emitter is also a good absorber). The Stefan-Boltzmann constant σ also cancels, and we just get T_E⁴ = 2T_A⁴. If we take the forth root of each side of the equation, we find that T_A = 0.84 T_E. The atmosphere is somewhat cooler than the earth.

Next, let’s reanalyze the power into and out of the earth when surrounded by an atmosphere. The sunlight power per unit area impinging on earth is still (1 – a)S/4. The radiation emitted by the earth is still σT_E⁴. However, the thermal radiation produced by the atmosphere that is aimed inward toward the earth is all absorbed by the earth (since the emissivity of the earth is one, e_E = 1), so this provides another factor of e_AσT_A⁴. Balancing power in and out gives

(1 – a)S/4 + e_AσT_A⁴ = σT_E⁴ .

Notice that if e_A were zero, this would be the same relationship as we found when there was no atmosphere: (1 – a)S/4 = σT_E⁴. The atmosphere provides additional heating, warming the earth.

We found earlier that T_E⁴ = 2T_A⁴. If we rewrite this as T_A⁴ = T_E⁴/2 and plug that into our energy balance equation, we get

(1 – a)S/4 + e_AσT_E⁴/2 = σT_E⁴ .

With a bit of algebra, we find

(1 – a)S/4 = σT_E⁴ (1 – e_A/2) .

Solving for the earth’s temperature gives

T_E = ∜[(1 – a)S/4σ] ∜[1/(1 – e_A/2) ] .

If e_A were zero, this would be exactly the relationship we had for no atmosphere. The fraction of energy absorbed by the atmosphere is not zero, however, but is approximately e_A = 0.8. The atmosphere provides a dimensionless correction factor of ∜[1/(1 – e_A/2)]. The temperature we found previously, 254 K, is corrected by this factor, 1.136. We get T_E = 288.5 K = 15.5 °C. This is approximately the average temperature of the earth. Our atmosphere raised the earth’s temperature from –19°C to +15.5°C, a change of 34.5°C.

Climate Change

To understand climate change, we need to look more deeply into the meaning of the factor e_A, the fraction of energy absorbed by the atmosphere. The main constituents of the atmosphere—oxygen and nitrogen—are transparent to both visible and thermal radiation, so they don’t contribute to e_A. Thermal energy is primarily absorbed by greenhouse gases. Examples of such gases are water vapor, carbon dioxide, and methane. Methane is an excellent absorber of thermal radiation, but its concentration in the atmosphere is low. Water vapor is a good absorber, but water vapor is in equilibrium with liquid water, so it isn’t changing much. Carbon dioxide is a good absorber, has a relatively high concentration, and is being produced by burning fossil fuels, so a lot of our discussion about climate change focuses on carbon dioxide.

The key to understanding climate change is that greenhouse gasses like carbon dioxide affect the fraction of energy absorbed, e_A. Suppose an increase in the carbon dioxide concentration in the atmosphere increased e_A slightly, from 0.80 to 0.81. The correction factor  ∜(1/(1 – e_A/2) ) would increase from 1.136 to 1.139, changing the temperature from 288.5 K to 289.3 K, implying an increase in temperature of 0.8 K. Because changes in temperature are the same if expressed in kelvin or Celsius, this is a 0.8°C rise. A small change in e_A causes a significant change in the earth’s temperature. The more carbon dioxide in the atmosphere, the greater the temperature rise: Global warming.

Feedback

We have assumed the earth’s albedo, a, is a constant, but that is not strictly true. The albedo depends on how much snow and ice cover the earth. More snow and ice means more reflection, a larger albedo, a smaller amount of sunlight absorbed by the earth, and a lower temperature. But a lower temperature means more snow and ice. We have a viscous cycle: more snow and ice leads to a lower temperature which leads to more snow and ice, which leads to an even lower temperature, and so on. Intermediate Physics for Medicine and Biology dedicates an entire chapter to feedback, but it focuses mainly on negative feedback that tends to maintain a system in equilibrium. A viscous cycle is an example of positive feedback, which can lead to explosive change. An example from biology is the upstroke of a nerve action potential: an increase in the electrical voltage inside a nerve cell leads to an opening of sodium channels in the cell membrane, which lets positively charged sodium ions enter the cell, which causes the voltage inside the cell to increase even more. The earth’s climate has many such feedback loops. They are one of the reasons why climate modeling is so complicated.

Conclusion

Today I presented a simple description of the earth’s temperature and the impact of climate change. Many things were left out of this toy model. I ignored differences in temperature over the earth’s surface and within the atmosphere. I neglected ocean currents and the jet stream that move heat around the globe. I did not account for seasonal variations, or for other greenhouse gasses such as methane and water vapor, or how the amount of water vapor changes with temperature, or how clouds affect the albedo, and a myriad of other factors. Climate modeling is a complex subject. But toy models like I presented today provide insight into the underlying physical mechanisms. For that reason, they are crucial for understanding complex phenomena such as climate change.

A Continuum Model for Volume and Solute Transport in a Pore

As Gene Surdutovich and I prepare the 6th edition of Intermediate Physics for Medicine and Biology, we have to make many difficult decisions. We want to streamline the book, making it shorter and more focused on key concepts, with fewer digressions. Yet, what one instructor may view as “fat” another may consider part of the “meat.” One of these tough choices involves Section 5.9 (A Continuum Model for Volume and Solute Transport in a Pore).

Neither Gene nor I cover the rather long Sec. 5.9 when we teach our Biological Physics class; there just isn’t enough time. So, at the moment this section has been axed from the 6th edition. It now lies abandoned on the cutting room floor. (But, using LaTex’s “comment” feature we could reinstate it in a moment; there’s always hope.) Russ Hobbie would probably object, because I know he was fond of that material. Today, I want to revisit that section once more, for old times sake.

The section develops a model of solute flow through a pores in a membrane. One key parameter it derives is the “reflection coefficient,” σ, which accounts for the size of the solute particle. If the solute radius, a, is small compared to the pore radius, R_p, then solute can easily pass through and almost none is “reflected” or excluded from passing through the pore. In that case, the reflection coefficient goes to zero. If the solute radius is larger than the pore radius, the solute can’t pass through (it’s too big!); it’s completely blocked and the reflection coefficient is one. The transition from σ = 0 to σ = 1 for medium-sized solute particles depends on the pore model.

The fifth edition of IPMB presents two models to calculate how the reflection coefficient varies with solute radius. The figure below summarizes them. It is similar to Fig. 5.15 in IPMB, but is drawn with Mathematica as many of the figures in the 6th edition will be. 

The blue curve shows σ as a function of ξ = a/R_p, and represents the “steric factor” 2ξ − ξ². It arises from a model that assumes there is plug flow of solvent (usually water) through the pore; the flow velocity does not depend on position. The maize curve shows a more complex model that accounts for Poiseuille flow in the pore (no flow at the pore edge and a parabolic flow distribution that peaks in the pore center), and gives the reflection coefficient as 4ξ² – 4ξ³ + ξ⁴. (Is it a coincidence that I use the University of Michigan’s school colors, blue and maize, for the two curves? Actually, it is.) Both vary between zero and one.

You can consult the textbook for the mathematical derivations of these functions. Today, I want to see if we can understand them qualitatively. For plug flow, reflection occurs if the solute is within one particle radius of the pore edge. In that case, the number of particles that reflect grows linearly with particle radius. The steric factor 2ξ − ξ² has this behavior. For Poiseuille flow, the size of the particle relative to the pore radius similarly plays a role. However, the flow is zero near the pore wall. Therefore, tiny particles adjacent to the edge did not contribute much to the flow anyway, so making them slightly larger does not make much difference. The reflection coefficient grows quadratically near ξ = 0, because as the particle radius increases you have more particles that would be blocked by the pore edge, and because the larger size of the particle means that it experiences a greater flow of solvent as you move radially in from the pore edge. So, the relative behavior of the two curves for small radius makes sense. In fact, for small values of ξ the two functions are quite different. At ξ = 0.1, the blue curve is over five times larger than the maize curve.

I find explaining what is happening for ξ approximately equal to one is more difficult. For plug flow, when the solute particle is just slightly smaller than the pore radius, it barely fits. But for Poiseuille flow, the particle not only barely fits, but it blocks all the fast flow near the pore center and you only get a contribution from the slow flow near the edge. This causes the maize curve to be more sensitive to what is happening near ξ = 1 than is the blue curve. I don’t find this explanation as intuitively obvious as the one in the previous paragraph, but it highlights an approximation that becomes important near ξ = 1. The model does not account for adjustment of the flow of solvent when the solute particle is relatively large are disrupts the flow. This can’t really be true. If a particle almost plugged a pore, it must affect the flow distribution. I suspect that the Poiseuille model is most useful for small values of ξ, but the behavior at large ξ (near one) should be taken with a grain of salt.

If find that it’s useful to force yourself (or your student) to provide physical interpretations of mathematical expressions, even when they’re not so obvious. Remember, the goal of doing these analytical toy models is to gain insight.

For those of you who might be disappointed to see Section 5.9 go, my advice is don’t toss out your 5^th edition when you buy the 6^th (and I’m assuming all of my dear readers will indeed buy the 6^th edition). Stash the 5^th edition away in your auxiliary bookshelf (or donate it to your school library), and pull it out if you really want a good continuum model for volume and solute transport in a pore.

Extracellular Magnetic Measurements to Determine the Transmembrane Action Potential and the Membrane Conduction Current in a Single Giant Axon

Forty years ago today I was attending my first scientific meeting: The Society for Neuroscience 14^th Annual Meeting, held in Anaheim, California (October 10–15, 1984). As a 24-year-old graduate student in the Department of Physics at Vanderbilt University, I presented a poster based on the abstract shown below: “Extracellular Magnetic Measurements to Determine the Transmembrane Action Potential and the Membrane Conduction Current in a Single Giant Axon.”

I can’t remember much about the meeting. I’m sure I flew to California from Nashville, Tennessee, but I can’t recall if my PhD advisor John Wikswo went with me (his name is not listed on any meeting abstract except the one we presented). I believe the meeting was held at the Anaheim Convention Center. I remember walking along the sidewalk outside of Disneyland, but I didn’t go in (I had visited there with my parents as a child).

Neuroscience Society meetings are huge. This one had over 300 sessions and more than 4000 abstracts submitted. In the Oct. 11, 1984 entry in my research notebook, I wrote “My poster session went OK. Several people were quite enthusiastic.” I took notes from talks I listened to, including James Hudspeth discussing hearing, a Presidential Symposium by Gerald Fischbach, and a talk about synaptic biology and learning by Eric Kandel. I was there when Theodore Bullock and Susumu Hagiwara were awarded the Ralph W. Gerard Prize in Neuroscience.

The research Wikswo and I reported in our abstract was eventually published in my first two peer-reviewed journal articles:

Barach, J. P., B. J. Roth and J. P. Wikswo, Jr., 1985, Magnetic Measurements of Action Currents in a Single Nerve Axon: A Core-Conductor Model. IEEE Transactions on Biomedical Engineering, Volume 32, Pages 136-140.

Roth, B. J. and J. P. Wikswo, Jr., 1985, The Magnetic Field of a Single Axon: A Comparison of Theory and Experiment. Biophysical Journal, Volume 48, Pages 93-109.

Both are cited in Chapter 8 of Intermediate Physics for Medicine and Biology.

This neuroscience abstract was not my first publication. I was listed as a coauthor on an abstract to the 1983 March Meeting of the American Physical Society, based on some research I helped with as an undergraduate physics major at the University of Kansas. But I didn’t attend that meeting. In my CV, I have only one publication listed for 1983 and one again in 1984. Then in 1985, they started coming fast and furious. 

Four decades is a long time, but it seems like yesterday.

The Difference between Traditional Magnetic Stimulation and Microcoil Stimulation: Threshold and the Electric Field Gradient

In Chapter 7 of Intermediate Physics for Medicine and Biology, Russ Hobbie and I discuss electrical stimulation of nerves. In particular, we describe how neural excitation depends on the duration of the stimulus pulse, leading to the strength-duration curve.

The strength-duration curve for current was first described by Lapicque (1909) as

where i is the current required for stimulation, i_R is the rheobase [the minimum current required for a long stimulus pulse], t is the duration of the pulse, and t_c is chronaxie, the duration of the pulse that requires twice the rheobase current.

An axon is difficult to excite using a brief pulse, and you have to apply a strong current. This behavior arises because the axon has its own characteristic time, τ (about 1 ms), which is basically the resistance-capacitance (RC) time constant of the cell membrane. If the stimulus duration is much shorter than this time constant, the stimulus strength must increase.

A nerve axon not only has a time constant τ, but also a space constant λ. Is there a similar spatial behavior when exciting a nerve? This is the question my graduate student Mohammed Alzahrani and I addressed in our recent article “The Difference between Traditional Magnetic Stimulation and Microcoil Stimulation: Threshold and the Electric Field Gradient” (Applied Sciences, Volume 14, Article 8349, 2024). The question becomes important during magnetic stimulation with a microcoil. Magnetic stimulation occurs when a pulse of current is passed through a coil held near the head. The changing magnetic field induces an electric field in the brain, and this electric field excites neurons. Recently, researchers have proposed performing magnetic stimulation using tiny “microcoils” that would be implanted in the brain. (Will such microcoils really work? That’s a long story, see here and here.) If the coil is only 100 microns in size, the induced electric field distribution will be quite localized. In fact, it may exist over a distance that’s short compared to the typical space constant of a nerve axon (about 1 mm). Mohammed and I calculated the response of a nerve to the electric field from a microcoil, and found that for a localized field the stimulus strength required for excitation is large.

Figure 6 of our article, reproduced below, plots the gradient of the induced electric field dE_x/dx (which, in this case, is the stimulus strength) versus the parameter b (which characterizes the spatial width of the electric field distribution). Note that unlike the plot of the strength-duration curve above, Fig. 6 is a log-log plot. 

Figure 6 from Alzahrani and Roth, Appl. Sci., 14:8349, 2024

We wrote

Our strength-spatial extent curve in Figure 6 for magnetic stimulation is analogous to the strength-duration curve for electrical stimulation if we replace the stimulus duration [t] by the spatial extent of the stimulus b and the time constant τ by the [space] constant λ. Our results in Figure 6 have a “spatial rheobase” dE_x/dx value (1853 mV/cm²) for large values of spatial extent b. At small values of b, the value of dE_x/dx rises. If we wanted to define a “spatial chronaxie” (the value of b for which the threshold value of dE_x/dx rises by a factor of two), it would be about half a centimeter.

To learn more about this effect you can read our paper, which was published open access so its available free to everyone. Some researchers have used a value of dE_x/dx found when stimulating with a large coil held outside the head to estimate the threshold stimulus strength using a microcoil. We ended the paper with this warning:

In conclusion, our results show that even in the case of long, straight nerve fibers, you should not use a threshold value of dE_x/dx in a microcoil experiment that was obtained from a traditional magnetic stimulation experiment with a large coil. The threshold value must be scaled to account for the spatial extent of the dE_x/dx distribution. Magnetic stimulation is inherently more difficult for microcoils than for traditional large coils, and for the same reason, electrical stimulation is more difficult for short-duration stimulus pulses than for long-duration pulses. The nerve axon has its own time and space constants, and if the pulse duration or the extent of the dE_x/dx distribution is smaller than these constants, the threshold stimulation will rise. For microcoil stimulation, the increase can be dramatic.