How to Solve a Nonlinear System

In a linear system the expressions, yp, for rates of change are linear functions of the dependent variables, y. More complicated functions of y do not appear, not even products of the dependent variables like y1 * y2. But most theoretical problems in Earth system science involve nonlinearities. For example, the rate at which a chemical reaction consumes species 1 may be proportional to the product of the concentrations of species 1 and the species 2 with which it is reacting, y1 * y2. In this chapter I shall describe and solve a simple nonlinear system involving the reaction between dissolved oxygen and organic carbon in the deep sea. I shall show how the nonlinear system can be represented by a linear system, provided that changes in the dependent variables are made in sufficiently small increments. Such increments are kept small by stepping forward in time with small steps. The time step can be adjusted automatically during the calculation so as to keep the increments small enough but no smaller than necessary. Steps that are too large cause errors or even instability, and steps that are too small waste time. The representation of a nonlinear system by its linear equivalent (for small increments) calls on algebraic manipulations that can be tedious and a prolific source of mistakes in complicated systems. This algebra can be avoided, however, by letting the computer perform the equivalent manipulations numerically. I shall demonstrate how to do this, finishing the chapter with a program that can solve coupled nonlinear systems directly from the equations for the rates of change of the dependent variables, automatically adjusting the time step to small values when the rates of change are large and to large values when the rates of change are small. Figure 4-1 shows a simple model of the processes that control the oxygen balance of the deep sea. Dissolved oxygen in surface seawater equilibrates with the atmosphere. Downwelling currents carry dissolved oxygen into the deep sea, where most of it reacts (metabolically) with organic matter carried into the deep sea in the form of settling particles of biological origin.

Climate: A Chain of Identical Reservoirs

One class of important problems involves diffusion in a single spatial dimension, for example, height profiles of reactive constituents in a turbulently mixing atmosphere, profiles of concentration as a function of depth in the ocean or other body of water, diffusion and diagenesis within sediments, and calculation of temperatures as a function of depth or position in a variety of media. The one-dimensional diffusion problem typically yields a chain of interacting reservoirs that exchange the species of interest only with the immediately adjacent reservoirs. In the mathematical formulation of the problem, each differential equation is coupled only to adjacent differential equations and not to more distant ones. Substantial economies of computation can therefore be achieved, making it possible to deal with a larger number of reservoirs and corresponding differential equations. In this chapter I shall explain how to solve a one-dimensional diffusion problem efficiently, performing only the necessary calculations. The example I shall use is the calculation of the zonally averaged temperature of the surface of the Earth (that is, the temperature averaged over all longitudes as a function of latitude). I first present an energy balance climate model that calculates zonally averaged temperatures as a function of latitude in terms of the absorption of solar energy, which is a function of latitude, the emission of long-wave planetary radiation to space, which is a function of temperature, and the transport of heat from one latitude to another. This heat transport is represented as a diffusive process, dependent on the temperature gradient or the difference between temperatures in adjacent latitude bands. I use the energy balance climate model first to calculate annual average temperature as a function of latitude, comparing the calculated results with observed values and tuning the simulation by adjusting the diffusion parameter that describes the transport of energy between latitudes. I then show that most of the elements of the sleq array for this problem are zero. Nonzero elements are present only on the diagonal and immediately adjacent to the diagonal. The array has this property because each differential equation for temperature in a latitude band is coupled only to temperatures in the adjacent latitude bands.

How to Deal with Several Reservoirs

The previous chapter showed how the reverse Euler method can be used to solve numerically an ordinary first-order linear differential equation. Most problems in geochemical dynamics involve systems of coupled equations describing related properties of the environment in a number of different reservoirs. In this chapter I shall show how such coupled systems may be treated. I consider first a steady-state situation that yields a system of coupled linear algebraic equations. Such a system can readily be solved by a method called Gaussian elimination and back substitution. I shall present a subroutine, GAUSS, that implements this method. The more interesting problems tend to be neither steady state nor linear, and the reverse Euler method can be applied to coupled systems of ordinary differential equations. As it happens, the application requires solving a system of linear algebraic equations, and so subroutine GAUSS can be put to work at once to solve a linear system that evolves in time. The solution of nonlinear systems will be taken up in the next chapter. Most simulations of environmental change involve several interacting reservoirs. In this chapter I shall explain how to apply the numerical scheme described in the previous chapter to a system of coupled equations. Figure 3-1, adapted from Broecker and Peng (1982, p. 382), is an example of a coupled system. The figure presents a simple description of the general circulation of the ocean, showing the exchange of water in Sverdrups (1 Sverdrup = 106 m3/sec) among five oceanic reservoirs and also the addition of river water to the surface reservoirs and the removal of an equal volume of water by evaporation. The problem is to calculate the steady-state concentration of dissolved phosphate in the five oceanic reservoirs, assuming that 95 percent of all the phosphate carried into each surface reservoir is consumed by plankton and carried downward in particulate form into the underlying deep reservoir. The remaining 5 percent of the incoming phosphate is carried out of the surface reservoir still in solution.

Why Simulate?

Our world is a product of complex interactions among atmosphere, ocean, rocks, and life that Earth system science seeks to understand. Earth system science deals with such properties of the environment as composition and climate and populations and the ways in which they affect one another. It also concerns how these interactions caused environmental properties to change in the past and how they may change in the future. The Earth system can be studied quantitatively because its properties can be represented by numbers. At present, however, most of the numbers in Earth system science are observational rather than theoretical, and so our description of the Earth system's objective properties is much more complete than our quantitative understanding of how the system works. Quantitative theoretical understanding grows out of a simulation of the system or parts of the system and numerical experimentation with simulated systems. Simulation experiments can answer questions like What is the effect of this feature? or What would happen in that situation? Simulation also gives meaning to observations by showing how they may be related. As an illustration, consider that area of Earth system science known as global change. There is now an unambiguous observational record of global change in many important areas of the environment. For elements of climate and atmospheric composition this record is based on direct measurement over periods of a decade to a century. For other environmental variables, particularly those related to the composition of the ocean, the record of change consists of measurements of isotopic or trace-element composition of sediments deposited over millions of years. This evidence of global change is profoundly affecting our view of what the future holds in store for us and what options exist. It should also influence our understanding of how the interaction of biota and environment has changed the course of Earth history. But despite the importance of global change to our prospects for the future and our understanding of the past, the mechanisms of change are little understood. There are many speculative suggestions about the causes of change but few quantitative and convincing tests of these suggestions.

Interacting Species in Identical Reservoirs

In Chapter 7 I showed how much computational effort could be avoided in a system consisting of a chain of identical equations each coupled just to its neighboring equations. Such systems arise in linear diffusion and heat conduction problems. It is possible to save computational effort because the sleq array that describes the system of simultaneous linear algebraic equations that must be solved has elements different from zero on and immediately adjacent to the diagonal only. This general approach works also for one-dimensional diffusion problems involving several interacting species. In such a system the concentration of a particular species in a particular reservoir is coupled to the concentrations of other species in the same reservoir by reactions between species and is coupled also to adjacent reservoirs by transport between reservoirs. If the differential equations that describe such a system are arranged in appropriate order, with the equations for each species in a given reservoir followed by the equations for each species in the next reservoir and so on, the sleq array still will have elements different from zero close to the diagonal only. All the nonzero elements lie no farther from the diagonal than the number of species. More distant elements are zero. Again, much computation can be eliminated by taking advantage of this pattern. I will show how to solve such a system in this chapter, introducing two new solution subroutines, GAUSSND and SLOPERND, to replace GAUSSD and SLOPERD. I shall apply the new method of solution to a problem of early diagenesis in carbonate sediments. I calculate the properties of the pore fluid in the sediment as a function of depth and time. The different reservoirs are successive layers of sediment at increasing depth. The fluid's composition is affected by diffusion between sedimentary layers and between the top layer and the overlying seawater, the oxidation of organic carbon, and the dissolution or precipitation of calcium carbonate. Because I assume that the rate of oxidation of organic carbon decreases exponentially with increasing depth, there must be more chemical activity at shallow depths in the sediment than at great depths.

How to Calculate Isotope Ratios

The calculation of isotope ratios requires special consideration because isotope ratios, unlike matter or energy, are not conserved. In this chapter I shall show how extra terms arise in the equations for the rates of change of isotope ratios. The equations developed here are quite general and can be applied to most of the isotope systems used in geochemistry. As an example of the application of these new equations, I shall demonstrate a simulation of the carbon isotopic composition of ocean and atmosphere and then use this simulation to examine the influence on carbon isotopes of the combustion of fossil fuels. As an alternative application I shall simulate the carbon isotopic composition of the water in an evaporating lagoon and show how the composition and other properties of this water might be affected by seasonal changes in evaporation rate, water temperature, and biological productivity. Equations for the rates of change of individual isotopes in a reservoir are not essentially different from the equations for the rates of change of chemical species. Isotopic abundances, however, are generally expressed as ratios of one isotope to another and, moreover, not just as the ratio but also as the departure of the ratio from a standard. This circumstance introduces some algebra into the derivation of an isotopic conservation equation. It is convenient to pursue this algebra just once, as I shall in this section, after which all isotope simulations can be formulated in the same way. I shall use the carbon isotopes to illustrate this derivation, but the same approach can be used for the isotopes of other elements, such as sulfur, oxygen, nitrogen, hydrogen, or strontium. The most abundant isotope of carbon has a mass of 12 atomic mass units, 12C. A less abundant stable isotope is 13C. And much less abundant is the radioactive isotope 14C, also called radiocarbon. It is convenient to express the abundances of these rare isotopes in terms of ratios of the number of atoms of the rare isotope in a sample to the number of atoms of the abundant isotope. We call this ratio r, generally a very small number.

The Carbon System and Several Useful Procedures

The routines developed in previous chapters can be used to simulate a variety of interesting systems in geochemical dynamics and global change. In addition to these routines are devices and procedures that can make easier the process of developing and debugging a simulation. I shall present several such procedures in this chapter, including the management of input and output files, particularly files of starting values, the definition of mnemonic names for variables, a graphic subroutine that provides a runtime view of the progress of a calculation, and the specification of complicated histories by means of a table. These are procedures that I find helpful, but because working with a small computer is a personal matter, you may not find them helpful. By all means, develop your own procedures or modify mine. As an application of these computational helpers I shall also introduce the carbon system and the equilibrium relationships among the species of carbon dissolved in natural waters. Carbon dioxide plays a key role in climate, in biological processes, in weathering reactions, and in marine chemistry. I shall next describe how the partial pressure of this gas in the atmosphere may be calculated. Because there is a rapid exchange of carbon dioxide between ocean and atmosphere, we must consider the fate of dissolved carbon. Carbon dissolved in seawater takes part in fast chemical reactions involving the species dissolved carbon dioxide H2CO3, bicarbonate ions HCO-3, and carbonate ions CO=3. The concentrations of these species are governed by equilibrium relationships (Broecker and Peng, 1982).

How to Calculate a Compositional History

The most interesting theoretical problems in Earth system science cannot be solved by analytical methods; their solutions cannot be expressed as algebraic expressions; and so numerical solutions are needed. In this chapter I shall introduce a method of numerical solution that can be applied to a wide range of simulations and yet is easy to use. In later chapters I shall elaborate and apply this method to a variety of situations. All numerical solutions of differential equations involve some degree of approximation. Derivatives—properly defined in terms of infinitesimally small increments—are approximated by finite differences: dy/dx is approximated by dely/delx, whose accuracy increases as the finite difference, delx, is reduced. A large value of delx may even cause numerical instability, yielding a numerical solution that is altogether different from the true solution of the original system of differential equations. But a small value of delx generally implies that many steps must be taken to evolve the solution from a starting value of x to a finishing value. A numerical solution, therefore, requires a trade-off between computational speed and accuracy. We seek an efficient and stable method of calculation that provides accuracy appropriate to our knowledge of the physical system being simulated. The problem described in this chapter can be easily solved analytically, and the analytical solution serves as a check on the accuracy of the numerical solutions. As a simple example of a global geochemical simulation, consider the exchange of carbon dioxide between the ocean and the atmosphere. The atmosphere contains 5.6 × 1016 moles of carbon dioxide (cf. Walker, 1977), a quantity that I assume to be in equilibrium with the ocean. In this illustration I assume that the oceanic reservoir is very large and therefore does not change with time. According to Broecker and Peng (1982, p. 680) the annual exchange of carbon dioxide between ocean and atmosphere is 6.5 × 1015 moles. The rate of transfer from atmosphere to ocean is proportional to the amount in the atmosphere; the flow from ocean to atmosphere is constant. Figure 2-1 summarizes this system. The amount of carbon dioxide in the atmosphere is proportional to the partial pressure.