Introducing key theoretical and data analysis tools in computational physics via Earth's temperature and climate

Teaching computation prepares our undergraduate students not only for graduate school but also for private sector careers, many of which will be in fields related to computation, such as programming and data science.^1,2 A perennial question facing teachers is how much to make our students do “by hand.” In computation, we have long made extensive use of programs written by others. This borrowing is efficient, but can lead to black boxes, where users do not know how the program actually works. The rise of machine learning, including artificial intelligence, has dramatically increased this black-box effect.

As physics teachers, we care about results, but we also care about derivations, because we care about understanding. We provide and ask for derivations in our analytical courses, and in traditional computation, we ask our students to write common algorithms prior to letting them use existing packages. There is somewhat less standardization when it comes to introducing students to data analysis, but it can be useful to have a broadly applicable and standardized approach here as well. Data and models about local temperature and climate change are well suited to this task, because they are topics about which students already have considerable intuition, while still providing enough complexity to thoroughly engage students.

The pair of activities I describe in the following, one in theory and one in data analysis, are linked only thematically and can be used separately, including with students at a wide variety of levels. My own target audience has been a broad group of undergraduate students who have all completed first-year physics but have little else in common. Many are physics majors, but others are majoring in math, chemistry, engineering, or other fields. About half of them have some prior programming experience, but the remaining half do not, and thus as a prerequisite to these activities, I spend a few weeks introducing them to Python. Any very high-level language suitable for scientific programming can be used, but Python stands out as one that is widely used in both scientific and corporate environments.

Adaptations can be made for students at different levels, as I discuss in Sec. IV. These activities can be used together or separately. Jupyter Notebooks implementing these activities in Python, including full solutions, are provided as supplementary material. These notebooks can be run locally on computers with Python and Jupyter installed (such as with the Anaconda distribution), or on a cloud-based service such as Google Colab.

Physics is stated in differential equations. While some situations are exactly solvable, other situations are not, particularly as students progress to more advanced and constrained versions of physical laws, such as Schrödinger's equation or Einstein's equations. Thus, we often resort to numerical solutions.

However, the equations of physics are almost always second-order equations, from the examples above to Maxwell's equations, Lagrange's equations, and even Newton's second law. Computational approaches to second-order equations can require learning how to simultaneously solve two first-order equations and then conceptualizing one second-order equation as two first-order equations.³ While this is an important skill, it is a less-than-ideal instructional starting point.

However, there are situations that give rise to first-order equations, including temperature change resulting from an imbalance of energy flow. A salient example is climate change.

This exercise also reviews physics concepts and introduces the idea of a toy model. Importantly, for numerical simulations, it also introduces validation of programs by checking against known results, and the concept of a converged simulation.

I start with the standard equilibrium calculation of a planet's temperature. Such calculations are commonly covered when training climate scientists^4,5 and astronomers,⁶ but do not actually require any specialized knowledge. First-year physics understanding of radiative heat transfer, based on the Stefan–Boltzmann equation, is all that is needed.^7,8

I do discuss the dependence on emissivity (⁠ $ϵ$⁠) and transmission coefficient (here denoted by R, since T is taken by temperature), but since these factors only occur multiplied together in this simple model, I instruct the students to combine these degenerate variables into a single effective transmission coefficient $Reff$⁠. As a pedagogical matter, it is worth emphasizing that the emissivities at visible wavelengths (which is most of the incoming radiation) and in the IR (most of the outgoing radiation) are different, as are the atmospheric transmission coefficients. More advanced students can appreciate that we are taking averages over more complicated, frequency-dependent transmission coefficients,^8,9 just as the Stefan–Boltzmann law averages blackbody radiation over all frequencies.

The first program the students write applies this analytical model to calculate equilibrium temperature for given transmission coefficients. This allows them to confirm that the transmission coefficients I provide do give a realistic global average temperature. It also lets them play with a variety of other scenarios, including the (low!) temperature of an Earth with no atmosphere, or the (also low) temperature of a “snowball Earth” that's quite reflective due to very large-scale glaciation, such as the Marinoan glaciation period around 650 million years ago.¹⁰ 

When the students write the non-equilibrium program, which does not incorporate Eq. (2), the first test they run with the program is to verify that it returns a stable temperature if they start with the correct equilibrium temperature for the given transmission coefficients. This introduces the concept of validating code by checking the results in known special cases. (Questions and solutions are shown in the climate model supplementary file.)

The non-equilibrium situation is introduced as a toy model for climate change. I have students employ a step-function decrease of $Routeff$ at five years into the simulation, a very simple model of increased outward opacity due to greenhouse gasses being introduced to the atmosphere. This abrupt change in transmission may seem unphysical, and I do primarily choose it for its simplicity, but note that sudden forcing can result from the climate system leaving quasi-stable equilibrium.⁴ The purpose of having this change five years into the simulation is to force another check that the non-equilibrium code can reproduce the initial equilibrium situation.

Students at this level are often unfamiliar with the term “toy model,” so it's worth emphasizing the utility of a model that captures some effect of interest even if it is not intended to realistically match data to any degree of accuracy. This is a single-zone model, using one temperature to characterize the Earth, and thus assumes perfect heat transfer around the globe. While clearly not the planet we live on, the toy model gives an example of how temperature increases can be “locked in” even if all new emissions of greenhouse gasses immediately cease.¹¹ 

The precise values of the transmission coefficients and heat capacity do not matter in this toy model, but make sure whatever combination you provide results in a correct initial equilibrium temperature, and a long-term temperature rise that is significant but not absurd. Instructors who have more time can introduce physically motivated values for these parameters. The heat capacity calculation can be motivated by taking the volume of water in the Earth's oceans,¹² about $1.3×1018$ m³, and calculating its heat capacity of about $5.4×1024$ J/(kg °C). The value I use in the code in my supplementary files is an order of magnitude smaller than this, and a more involved discussion could delve into the approximations we are making in this single-zone model (which assumes perfect mixing, something that does not happen in our oceans).

Motivating the transmission coefficients would be a more significant digression, appropriate only for an advanced class.^9,13 This could be worthwhile if the class has already covered numerical integration, in which case the students could integrate over an empirical atmospheric transmission function. You will need full coverage from the near UV through the thermal IR. Since no single instrument can record these data, the best approach might be an empirically based radiative transfer model, such as MODTRAN6,¹⁴ or other codes making use of the absorption database HITRAN.¹⁵ Also note that you must weight the actual transmission function by the emission spectrum to calculate a single effective number.

I have the students calculate how long it takes to reach the new equilibrium, to within 0.1 K, and they see a temperature rise due to greenhouse gas emissions from decades earlier (Fig. 1). I also have the students vary their time step, $Δt$⁠, and check their results. At this early stage in their computational careers, they are tempted to choose time steps that seem like round numbers to them, not always distinguishing the large difference that can exist between $Δt=1$ (second) and $Δt=1$ (year).

The meaning of “convergence” is that the model produces the solution we care about to the precision we need. Remind students that when solving a differential equation, the solution is a function, not a number, and in our case, we are looking for the temperature function $T(t)$⁠. This definition of convergence will likely be novel to the students, who will be more familiar with the idea of convergence, meaning that a function approaches a specific number for limiting values of the independent variable.¹⁶ 

Since we do not know the true solution, we can check for convergence by decreasing the time step $Δt$⁠. If one sees a significant change in $T(t)$ when decreasing $Δt$⁠, that means the simulation is not converged. Using smaller time steps makes the simulation more accurate (until numerical noise starts to dominate), but decreasing the time step means increasing the time it takes to run the simulation. Ideally, then, we should use time steps that are not much smaller than the level needed for convergence.

In most of our classes, we try to solve things exactly, and it can be a challenge for students to let go of this mindset and aim for results that are “good enough.” It is also novel for the students to check their results not by comparing to the “right” answer, but rather by comparing two models when we do not know in advance that either is right (to within our tolerance).

Students find that this is a relatively gentle introduction to numerical differential equation solution, given the simplicity of the first-order differential equation and the sufficiency of the Euler approximation for solving the equation. The Euler approach, as implemented in Eq. (3), is only first-order accurate in the time step, and it is usually too inefficient to use in practice, but here it works fine. This is beneficial, as the Euler approach is straightforward to relate to the differential equation itself, Eq. (4), although this is not to be expected to be immediately obvious to students. From this basis in numerical approaches, instructors can lead into a discussion of higher-order Runge–Kutta techniques,³ as well as techniques to solve second-order ordinary and partial differential equations.

To tie this back to physics, it is worth closing by asking the students to explain why larger time steps always lead to faster temperature rises.

Too often the data we initially give students to analyze is simplified, sanitized, or even mock data. There is use to this, to initially illustrate some points with no added complication, but I think we need to quickly move past this into working with messy data. Students often get to encounter such data in research projects, but it is useful to prepare students for this research experience with an exercise that is standardized, but still has data that are messy and can also be made unique for each campus.

To this end, I chose the temperature datasets that National Oceanic and Atmospheric Administration (NOAA) provides through Climate Data Online¹⁷ (CDO). Hourly temperature readings are available for weather stations across the United States, which provides data for more complex model fitting. In particular, there are diurnal variations in temperature and annual variations in temperature, and on top of that there might be the subtle signal of long-term rising temperatures from climate change. In addition, the data are imperfect and may not be formatted as we would like, providing students with an example of the complications of real-world data.

First show students how to read in the data and turn it into two arrays of time (t) and temperature (T), as a physicist would expect. For time, I use the common astronomical standard of the modified Julian date (MJD), which is the floating-point number of days since midnight on July 17, 1858. Other formats, such as seconds since the Unix epoch, can work, although the CDO data can extend back well before 1970 for many stations. Note that the datasets provided through CDO are good, but imperfect. For example, in the dataset I have provided as an example, the year 1996 is entirely missing. Perhaps this weather station had a period of mourning after the end of the Calvin and Hobbes comic strip.¹⁸ 

In the supplementary material, I provide the questions I pose and full solutions in the temperature and climate data analysis notebook, as well as an example data file for one particular weather station. I start this data file in 1973. While this location has data back to the 1940s, it is increasingly sparse, and the weather station actually changed location during this period. There are nearly half a million data points, but I chose a plain text format to make it easily human readable. More advanced classes can introduce more efficient file formats.

However, even with the full seven free parameters, students get extremely disappointing results from curve fitting with Eq. (5). Their plots show extremely poor matches to the data, as seen in the dot–dash line of Fig. 2.

At this point, I draw their attention back to two exercises I had them do earlier. The first was using brute-force calculation to find the least-squares best fit to a collection of lab measurements of the local acceleration due to gravity (g). To accomplish this, they generated an array of trial g values and calculated the $χ2$ for each one. They then found the lowest $χ2$ and selected the g value associated with it as their best estimate. This may seem like an absurd approach, as the optimal result is just the arithmetic mean of the data, but this exercise has several benefits. First, it allowed them to easily interpret $χ2$ minimization visually, as you can make a 2D plot of $χ2$ vs. trial g and see the clear minimum. Such plots become difficult and then impossible as the number of parameters increases. Second, they saw that the least-squares approach, done properly, produces an expected result that they can find independently. Third, they saw that it does not quite produce the expected result, because they only had a finite set of trial g values. They can increase the accuracy by using more trial g values, but they saw the tradeoff between accuracy and the cost of more memory and higher runtime.

The second exercise is to use brute-force calculation to do a linear fit to some data, calculating $χ2$ over a 2D parameter space of slope and intercept. They chose 1-D arrays of slope values and of intercept values and generated a 2D array of paired slopes and intercepts. They calculated the $χ2$ for each pair of slopes and intercepts, and then, as with the first problem, they found the minimum $χ2$ and selected the associated slope and intercept as their best fit. Again this method has no practical application, given that the linear problem can be exactly solved with matrix inversion, but it showed them the computational side of problems that cannot be solved exactly. The slopes and intercepts they found were good, but not quite right, because the sampling that they could do was somewhat course-grained in a 2D space. Here, the cost of increasing accuracy is much larger than in the first problem of finding g, because the cost of the calculation goes like $N2$ rather than like N, where N is the number of points in the 1D arrays. They could still plot $χ2$ as a function of the model parameters, and they saw that while some regions were easily ruled out, there was a valley of fairly similar $χ2$ values. While the slope and intercept were not degenerate, there was an approximate degeneracy. With an upward sloping line, a decrease in intercept can approximately be compensated for by an increase in slope. Thus, this exercise showed them the increased difficulty of brute-force minimization with increased dimensionality and also showed that minimization can be tricky when multiple parameters are involved and can change in correlated ways.

This background, while not necessary, equips the students to understand what had gone wrong with scipy.curve_fit(). They have $105$– $106$ data points and a seven-dimensional parameter space. Brute-force $χ2$ minimization is essentially impossible, certainly at a fine enough level to get adequately precise results. We discuss the concept of a local minimum, and how easy it is for an algorithm to get stranded far from the optimal point in parameter space.

I have them read the help page for scipy.curve_fit(), and they notice that you can specify initial guesses for parameters. After some thought, they realize they can give extremely good initial guesses for $ω1$ and $ω2$⁠, and reasonable guesses for $A1$⁠, $A2$⁠, and C. It is possible to guess the phases as well, but it is not necessary; with the information from the periods, scipy.curve_fit() is able to find a good fit. Students will discover that the routine does not obey physics conventions. For example, it often returns negative amplitudes, which students may only recently have had drilled out of them in first-year physics.

The model plot with the true best-fit parameters still can look wrong to students. On annual scales, the amplitude looks too small, for example, as seen with the model fit in Fig. 2. Plots that zoom in on few-day regions can help explain what is happening. Some time periods show very little diurnal variation (Fig. 3), while others show much more (Fig. 4).

I have students use their temperature model to predict the temperature for our last class meeting. This is not entirely trivial, as they have to make sure their model (in UTC) is correctly shifted to match our time zone, and many will shift the wrong direction. However, this is a fun test, the like of which students may not have experienced in their classes before, where they are predicting something that is genuinely unknown but can be tested easily and viscerally quickly thereafter. The first time my students did this, their prediction was closer to reality than what the National Weather Service predicted 24 h in advance. This was obviously luck, given that our model knows nothing about current weather conditions, but it is a good example for students about how hard prediction is with complicated systems—even advanced computer models taking into account recent data can be inaccurate for complicated systems with many physical scales.

While I use an explicit two-sine model, I do mention Fourier analysis to my students and stress the utility of such an approach for processes where they do not have good guesses at the periods involved. The power spectrum thus obtained, shown in Fig. 5, confirms the two periods that we knew a priori were present, gives some information on their relative strengths, and confirms that no other periods of any significance occur in the data. My class has no physics or math prerequisites that would cover Fourier analysis, so, for me, this is just a demonstration, but other instructors could modify this based on their students' backgrounds.

The slope they get seems small until they appreciate the units. With the weather station I use, we find $α=(6.81±0.14)× 10−5 °$ C/day in the fit to the residuals, but this leads to a 1.1 °C rise from 1973 through 2018. (When directly fitting the original data, we find a slightly different number, but it is consistent: $α=(6.74±0.14)×10−5$ °C/day.)

I stress to the students that a single temperature station is not enough to establish a global trend like climate change, and that satellite observations are crucial for whole-Earth temperature analysis. Nonetheless, it is nice that in this sample dataset, the rise is nearly exactly equal to the global mean temperature rise over this period.

These activities are best completed as a combination of in-class work and homework, but the exact proportion of each is flexible and can be adapted to the circumstances of particular courses.

The sample code provided has a specific amount of scaffolding that is appropriate for students who have completed first-year physics but not necessarily anything beyond that and have received a brief introduction to programming, plotting, and statistics (perhaps in the current class, or perhaps earlier). However, the amount of scaffolding is very flexible. Remove some of the supports, and the problems become more challenging—indeed, the data analysis one would be appropriate to challenge graduate students if it had very little scaffolding.

Scaffolding can also be increased to allow use with lower-level students. For example, an instructor could use the temperature residuals after fitting out annual and diurnal variation to discuss climate change with a class that has only seen linear regression. When asking students to write a function, an instructor who wants to provide an intermediate level of scaffolding can provide the docstring, specifying the input and output, and allow the students to create the requisite code.

The activities discussed here are also ripe targets for extending into further projects. The simplest extension of the climate model is to assume that the transmission coefficient changes slowly over time, to mimic more realistic increases of greenhouse gasses in the atmosphere. The effective transmission in (which includes the Earth's albedo) could be changed as a function of temperature, for example modeling the transition of sea ice to open sea with its lower albedo. Implementing models with multiple zones to the Earth, with different temperatures and albedos, would be an interesting but more challenging extension. Of course there are many other extensions possible along the road to a professional climate model, such as taking into account some of the effects of convection, and the atmosphere's nonlinear and wavelength-dependent opacity response to additional greenhouse gasses,¹³ accounting for feedback,⁷ creating a 1D atmospheric model,⁹ and so on.⁵ 

The temperature data have many extensions. For example, in my class, students are introduced to confidence intervals. Thus, I have the students calculate 95% confidence intervals on the average temperature at a specific day and time using both Gaussian statistics and a nonparametric bootstrap. (Daylight saving time can be nontrivial here, as some dates have been on standard time some years and daylight saving time other years.) They can also find normal temperature ranges that could be used in calculations of expected heating and cooling loads for buildings in the area. It would also be interesting to use data from many different weather stations to check for consistency on climate change, or compare the temperature stability of various regions.

The Fourier analysis, presented to my students just to let them know that such approaches exist, would form the beginning of the analysis for students familiar with such approaches. Do be careful here. Discrete Fourier transform algorithms tend to assume uniformly sampled data, and real weather station data are nearly but not quite uniform. In the example data I provide, 99.923% of the data are separated by 1 h, but you need 100%. I use spline interpolation, which my students have briefly seen before, to regularize my data.

Scientific data analysis is making increased use of machine learning, including AI, which presents opportunities but also pitfalls. Machine learning approaches create their own algorithms or statistical weights that were not coded by any person, replacing boxes that were black to most users (like scipy.curve_fit()) with boxes fairly black to all users. These tools are powerful, and necessary to make cutting-edge progress in science, but, as educators, we must work to lift the veil and give students at least a glimpse of what is going on behind the scenes. We must also make sure that students think critically about the results they get, from any routine or model, and build the skills to assess whether or not things make sense. This has become an issue of great importance in the era of generative AI, where probabilistic output means you can no longer say an algorithm is clearly correct or incorrect, but rather each output must be assessed independently.

Temperature data are well suited to introducing both theoretical techniques and data analysis because of its familiarity to the students. This reduces the cognitive load when it comes to the data and allows them to focus on the novel aspects of numerical simulation and data analysis.

Please click on this link to access the supplementary material online for a text data file containing temperature information for one weather station from 1973 to 2018. Two Jupyter Notebooks of text and Python code are also included. The first describes the climate model of Sec. II. The second goes through the data analysis and model fitting of Sec. III. Print readers can see the supplementary material at https://doi.org/10.60893/figshare.ajp.c.7523670.

The author has no conflicts to disclose.