by Dan Hughes and Tomas Milanovic
Further reflections on the application of the divergence theorem to the Earth’s climate system.
The subject of a recent post at CE was application of the divergence theorem to some aspects of Earth’s climate systems. Application of this theorem to an equation for conservation of the total energy of the atmosphere led to requirements for that energy to be constant. A corresponding post appeared at . . .aTTP not long after the CE post. The post at . . .aTTP was focused on disagreements with the post at CE.
These posts and associated lively discussions led to disagreements among the commenters on several mathematical aspects of the subject. The objective of the present post is to keep focus on some points of major disagreement: And hopefully, agreements.
Two of the issues for which there was massive disagreement were associated with the use of some well established concepts from mathematics and thermodynamics. In particular, known results and theorems from Partial Differential Equation (PDE) theory as well as from non-linear dynamics need to be strictly respected.
We directly address the points of major disagreement and have tightened our mathematical and thermodynamic arguments. The foundations of our agreements in mathematics and thermodynamics that were used in the previous post are further explained and expanded.
The divergence theorem states that the oriented flux of a vector field through a closed surface is equal to the volume integral of the divergence over the region inside the surface.
Physically this theorem is used for different conservation laws (mass, energy, charge) which take the general form of continuity equation:
where q is the density of some quantity, J is a vector field representing the flow of this quantity and S are sources/sinks of this quantity. In the case of energy which can neither be created nor destroyed we have:
Integrating this equation over some volume V with boundary B and applying the divergence theorem we obtain the well known law that “The variation of energy in a volume V is equal to the oriented flux through its boundary B“.
Obviously if the continuity equation is written for non-conserved quantities like radiation or heat, it is necessary to take into account all sources and sinks of this quantity before integrating in the volume V .
Application of the theorem to idealized sub-systems of Earth’s climate system leads to equations for the time-rate-of-change of the total energy content of the sub-systems. The resulting equations must identify physical phenomena and processes that can lead to changes in energy content, and the locations of these. Modeling and evaluations of these will determine if the net flow of energy is positive of negative (as per the theorem), and thus determine the time-rate-of-change in total energy.
Importantly, the equations also lead to the requirements that must be met in order for the total energy content to remain constant.
It is also to be noted that when this theorem is applied to energy, the results state only the obvious – energy is conserved – but gives no information about the dynamics of the system.
The divergence theorem is also known as Gauss’s theorem, or Ostrogradsky’s theorem, or the Gauss-Ostrogradsky theorem, and finally Stokes’s theorem. Establishing historical precedent is sometimes a tricky, and long-term, process.
Earth’s Climate System and Equilibrium
There can be no question that Earth’s climate system is a system that operates far from equilibrium. Here we use equilibrium in the universally established and accepted sense that the word is used in thermodynamics. In thermodynamics, equilibrium means absence of all gradients: Nothing can change over time, and nothing can change in space.
Flows of mass, momentum, energy, chemical, and biological processes are present in Earth’s climate system, and these operate over very wide ranges of time scales. These flows require the existence of gradients in driving potentials. Gradients in driving potentials absolutely ensure that the system is never in equilibrium.
In the universally accepted sense of equilibrium, the gradients will exist for all time and Earth’s climate system will always be far from equilibrium. The flows change on all time and space scales, and these changes will affect, to lesser and greater extents, the radiative energy transport properties and characteristics of Earth’s climate.
When thermodynamics is the subject of discussions the universally accepted nomenclature should form the basis of word meanings. In this regard, Earth’s climate system cannot and will never attain an equilibrium state. More generally it must be noted that any rotating sphere exchanging radiation with a point-like source is never in equilibrium because different points on its surface are at different temperatures and therefore temperature gradients exist at all time.
Given the impossibility of attaining an equilibrium state, other possibilities for the state of Earth’s climate systems may only be; (1) unsteady state, (2) stationary state, or (3) steady state.
The GCMs, and the evolving Earth System Models (ESMs), in which the coupled nature of the components of Earth’s climate system are mathematically modeled, are all set as an Initial-Boundary Value Problem (IBVP).
It should be very clear from the previous post that the focus was exclusively on GCMs. It is very clear that the post addressed GCM models that are based on the interactions of the many various components of Earth’s climate system. Other, special purpose, models were not addressed.
Initial-Boundary Value Problem (IBVP)
A Partial Differential Equation (PDE) problem is said to be well posed if
1. a solution to the problem exists 2. the solution is unique, and 3. the solution depends continuously on the problem data.
The solutions to PDEs are arbitrary functions. General and mathematically rigorous proofs of the general existence, regularity, and uniqueness requirements for solutions of PDEs exist only for some simple cases. In particular they do not exist for the Navier-Stokes equations. Those well-posedness requirements, however, do not prevent applications of the concept of an IBVP to the present issues. To set a well-posed problem the correct number and kind of boundary conditions are also required. As a zeroth-order requirement, we always count the number of dependent variables and of their derivatives compared to the number of equations available to ensure that the numbers are the same.
We all know that a first order PDE in a spatial variable x and a time variable t has a unique solution (if it admits one) if and only if the solution u satisfies a boundary constraint (u=f on some boundary B of an open set U) AND/OR an initial value constraint (u=g for t=0). Absence of either of these leads to an infinite number of arbitrary solution functions. Boundary constraintand boundary value are generally interchangeable.
A boundary condition should be a function of x and t if the value of the unknown function also depends on x and t. For instance if a value is prescribed at some location B then u(B,t) = g(x,t)where g(x,t) is the forcing function imposed at B and can vary with x on B.
A special remark about transformations of boundary conditions is here necessary. If one spatially averaged the boundary condition g(x,t) one would obtain h(t) as boundary condition because the averaging destroys the x dependence. It follows that if the average is used, the boundary condition now becomes u(B,t) = h(t). This is clearly a completely different boundary condition than u(B,t) = g(x,t). It follows obviously that solutions of the same PDE with the first boundary conditions will not be the same as solutions with the second boundary conditions. In this case the solutions, if they exist, would be different from the solutions of the original PDE and therefore invalid.
This situation frequently arises in a great many (almost all) numerical calculations, whether or not explicit averaging operators have been applied to the basic continuous equations. Averaged equations are the rule, and not the exception, for all complex problems. Discrete approximations are a form of implicit averaging. If a grid is refined to higher resolution, the boundary conditions can change as described above. For example, the topology, as seen by the discrete equations, of the boundary surface itself can be effectively changed. The location of a solid embedded in a flow field will also change unless the grid refinement is carried out in a manner that considers this effect. These lower-order numerical artifacts prevent the order of the numerical method from reaching the theoretical asymptotic range, and can in fact result in the application order of the method to be less than unity.
In order to simplify an otherwise unsolvable equation, it arises often in practice of fluid dynamics that a change of variable involving averages is done: u(x,t) = U(x,T) + u'(x,t) where U(x,T)is the temporal average of u(x,t) over some interval [t,t+T] with T constant and u'(x,t) is a “fluctuation” around the average postulated to be with a temporal average 0. The result of this transformation is that the equation is no more naturally closed and empirical ad hoc closures are necessary.
Another far reaching consequence is that this transformation only simplifies the original equation if it is postulated that the short time scales of u'(x,t) do not interact with the larger time scales of U(x,T).
This is typically not the case in spatio-temporally chaotic systems where the time scales interact and U must be written as U(x,T,u'(x,t)) which defeats the original purpose which was to “decouple” the short and long time scales. Of course a similar change of variable could be done with a spatial average instead of the temporal average but the consequences would be identical.
ILL-Posed Initial-Boundary Value Problem
It has turned out that calculated GCM results indicate that the equation systems that form the basis of GCMs have resulted in ILL-Posed Initial-Boundary Value Problems. The numerical solutions do not depend continuously on the initial data, so Requirement 3 above is not met. Additionally it is not known if a unique solution exists so that Requirement 2 can only be postulated.
Analyses of the continuous equation systems, discrete approximations to these, and the applied numerical solution methods to establish a theoretical basis for the observed calculated results is far beyond the state-of-the-art for spatio-temporal chaos.
The hyperbolic formulation of the hydrostatic balance model for the momentum equation in the vertical direction, the basis of the dynamical cores for the atmosphere and ocean in several GCMs, has been shown to have complex characteristics; the equations are not hyperbolic under the conditions for which they should be. The ramifications of this result include; (1) artificial, purely numerical, growth of all wavelengths of perturbations, and (2) questions related to boundary-condition specifications. Item (1) is especially non-physical relative to the high frequency, very short wavelengths that in the physical domain are responsible for dissipation that stabilizes fluid flows.
When discussion are focused on mathematical matters, the universally accepted nomenclature should form the basis of word meanings. In this regard, GCMs using the hydrostatic approximation represent ILL-Posed Initial-Boundary Value Problems.
The WG1AR5 Final Report includes Glossary III that defines the terms as used in the report. On chaos, the Glossary says:
“Chaotic A dynamical system such as the climate system, governed by nonlinear deterministic equations (see Nonlinearity), may exhibit erratic or chaotic behaviour in the sense that very small changes in the initial state of the system in time lead to large and apparently unpredictable changes in its temporal evolution. Such chaotic behaviour limits the predictability of the state of a nonlinear dynamical system at specific future times, although changes in its statistics may still be predictable given changes in the system parameters or boundary conditions.”
This definition applies to temporally chaotic systems where the dynamical parameters (pressure, temperature . . .) depend only on time and the dynamics are governed by a finite number of nonlinear ordinary differential equations (ODEs). A system such as the climate system is governed by nonlinear PDEs so that variations in space are as chaotic as variations in time. As the perturbations in spatio-temporal chaos propagate in space, which they cannot do in temporal chaos where space doesn’t exist, the dynamics of spatio-temporal chaos has little in common with the mere temporal chaos and is much more complex.
In particular the ergodic theorem which allows to connect temporal averages to the phase space averages is valid only for temporal chaos if it is ergodic (for example the famous Lorenz attractor) but has no equivalent for spatio-temporal chaos.
From there follows that there is no reason to expect a unique invariant probability distribution in spatio-temporal chaos which could be determined from boundary constraints alone. It is neither proven nor expected that spatial and temporal probability distributions are invariant and predictable on large time scales. A typical example of the variability of probability distributions on large time scales are the chaotic ice age and interglacial oscillations. These oscillations illustrate that it is not only impossible to predict their dynamics but it is also impossible to predict probabilities that this or that particular event takes place.
A current error for those not familiar with nonlinear dynamics is to exclusively assimilate the chaos in climate to the very short term fluid fluctuations like turbulence. While these short term fluctuations are indeed chaotic and lead to the unpredictability of the weather which is the base on which the climate is defined, the chaos in the climate is also due to nonlinear interactions at longer time scales and is not limited to fluids. Ice dynamics, clouds, biological phenomena and variations of astronomical parameters contribute to climatic spatio-temporal chaos yet they are neither represented by fluids nor by short time scales.
Here it has to be noted that boundary and initial conditions must be specified prior to solving any given PDEs. This is particularly difficult for the climate system. A boundary of a physical 3D system is necessarily a surface so giving a boundary condition means to prescribe the value of the unknown functions for every point and every time at this surface.
Yet these values are by definition unknown because they depend on the nature of the unknown function inside the boundary and therefore cannot be prescribed Ã priori. In the most favorable case when at least the spatial average of a parameter could be prescribed on the boundary it would still not allow to determine correctly the unknown functions because the spatial average operator destroys the spatial variability so that the boundary condition becomes u(x,t) = g(t) at the boundary B instead of the correct u(x,t) = g(x,t) at B. This has for consequence that the solutions computed with h(t) as boundary condition are invalid.
The phrase ‘boundary condition’ appears 47 times in the entire WG1AR5 report, is used in the Glossary but not defined there. The phrase ‘boundary value’ is used very infrequently, and also is not defined in the Glossary. The word ‘constraint’, and variations ‘constraints’ and ‘constrained’, are widely used in the WG1AR5 report. Constraint is not defined in the Glossary.
Boundary Value Problem (BVP)
The BVP characterization offered by Professor Steve Easterbrook is widely cited. Professor Easterbrook said:
For understanding climate, we no longer need to worry about the initial values, we have to worry about the boundary values. These are the conditions that constraint the climate over the long term: the amount of energy received from the sun, the amount of energy radiated back into space from the earth, the amount of energy absorbed or emitted from oceans and land surfaces, and so on.
The first of these is indeed a boundary condition. The second is not a boundary condition, and additionally is incorrect both mathematically and physically. Specifying the energy, or mass, leaving a system violates the fundamental principles of energy and mass conservation. The energy absorbed or emitted by material interior to the system are calculated by the model equations; they are not specified.
Computational Fluid Dynamics and GCMs
In many discussions of models and numerical methods and GCMs, the relationship between Computational Fluid Dynamics (CFD) and GCMs comes up. These, of course, are related as fluid motions are at the heart of the dynamics of the atmosphere and oceans.
Generally, CFD applications are very significantly advanced in critically important aspects compared to GCMs. Very significant research and development efforts have been applied to CFD models, methods, software, and applications over several decades. Fully Verified and Validated commercial-grade CFD software, ready for a variety of everyday engineering applications, is a mature market outside the framework of research and development.
GCMs on the other hand are, relatively, more or less under constant research and development. Earth’s climate system is vastly more inherently complex than all and any CFD applications because the latter deal uniquely with Navier-Stokes equations while the former must include, among others, chemical, biological and astrophysical changes as well as many critical phase change processes leading to clouds and ice cover dynamics.
Frequently within these discussions the fact that CFD is applied to various aspects of commercial aircraft is noted. Of course that’s correct. However, these CFD applications are not related in any way whatsoever to the safety of commercial aircraft flight. Moreover the commercial CFD application are certified only for precisely specified ranges of the dynamical parameters. Typically the aircraft behaviour in stall conditions which is fundamental for the safety is outside of the specified range of parameters.
The safety and flight worthiness of commercial aircraft are completely determined by the independent Federal Aviation Administration (FAA). Flight tests determine the safety and flight worthiness of commercial aircraft. Safety and flight worthiness of commercial aircraft are not functions of any CFD calculations. After all, aircraft had been determined to be flight worthy for decades before CFD entered the picture.
It is an obvious fallacy to say, “You’re safe flying in commercial aircraft”, and then conclude, “Therefore GCMs are correct.” CFD and GCMs work with very different systems but, as already mentioned, CFD for aircraft dynamics is used for a specified range of fundamental parameters that have been Validated in wind tunnels, while GCM are supposed to work for all ranges of parameters and cannot be tested in wind tunnels.
For the most part, almost all applications of CFD to aircraft focus on a very limited sub-set of fluid dynamics, mainly turbulence, using a limited number of Validated parameters that can be directly related to the physical phenomena and processes of interest, and the states of the fluid that will be encountered in the physical domain. Numerical methods are a second focus and are deeply investigated. Relative to spatial resolution, numerical grids for CFD applications per unit of total modeled volume contain more cells by several orders of magnitude than GCM applications.
Qout = Qin
Discussions at CE and . . .aTTP more or less centered around the concept that at some undefined future time, when averaged over the entire planet for some, also, undefined time period, the out-going radiative energy at the Top of the Atmosphere (ToA) will be equal to the in-coming radiative energy. This hypothesis seems to be the basis for the argument that the “final state”, also frequently characterized as an “equilibrium” state, is determined solely by “the boundary conditions”.
Note that we already know that an “equilibrium” state is never going to happen. Given that chaotic response is observed in GCM calculations, steady state, and stationary state cannot be applied, so unsteady state must describe GCM results.
Also note that to postulate Qout = Qin, even in average over a time T, it is necessary to postulate that no energy sinks or sources (phase changes, chemical reactions, storage) existed in the system during the time T. Additionally if the surface where Qin and Qout are defined is such that only radiation fluxes exist and one is interested only in the non-reflected fraction of the incoming radiation, it is necessary to postulate that the reflectivity stayed constant during the time T.
Last but not least we must strongly remind what has been already said above. Postulating “Qin = Qout” through some surface has for a trivial consequence that the total energy in some volume is constant (see Gauss theorem). This tautologically true statement has nothing to do with boundary conditions because it is true for any arbitrary volume, any arbitrary boundary, and any arbitrary value of the dynamical parameters for which “Qin=Qout” is postulated true.
In particular, this tautology gives no information whatsoever about the dynamics of the system where any arbitrary energy and momentum distribution may exist as long as the total energy is conserved and we may be reasonably sure that it will be indeed conserved everywhere and for all times.
The equality condition is used with some GCMs as the target while tuning the model by use of variations in the many parameterizations that are critically important to the degree of fidelity of calculated results with the physical domain.
Earth’s climate system is open with respect to energy; energy both enters and leaves the system. It is yet another universally accepted concept that the energy (and mass) leaving a thermodynamic system cannot ever be specified. The energy leaving a system cannot be set as a boundary condition, boundary value, or boundary constraint. This statement is correct no matter how far away is future time, or for how long the average is taken. The energy, or mass, leaving a system cannot be specified because that act in itself violates conservation of energy, or mass, and other foundational concepts of mathematics.
The argument seems to claim that energy conservation is the the only part of the system that matters for the dynamics. Initial conditions don’t matter, turbulence doesn’t matter, chaos doesn’t matter, ill-posedness doesn’t matter, parameterizations, and tweaking them, don’t matter, tuning doesn’t matter, discretization doesn’t matter, numerical solution methods don’t matter, application grid resolution doesn’t matter, inherent weaknesses of the continuous equations, especially the parameterizations, don’t matter, the known limitations and uncertainties of many of all these aspects don’t matter, among others.
An easy thought experiment which has never been done in practice to our knowledge would be to run a GCM with initial conditions specified as “Gulf stream exists at t=0” and “Gulf stream doesn’t exist at t=0” while the boundary conditions are identical for both cases. It is beyond any reasonable doubt that the computed solutions and their statistics (averages, standard deviations, and higher order momenta) would be extremely different, This experiment would represent a numerical demonstration that initial conditions matter on all time and space scales and for all parameters .
In contrast, considering the possibilities associated with the long list of aspects that don’t matter, it seems possible that these aspects have resulted in models, methods, software, and application procedures for which the initial conditions and the boundary conditions don’t matter.
This is the very foundation of the present post. So far as we are aware, the “boundary conditions” that solely determine the “final state” have never been identified. Neither have the future time at which this state is attained and the time-period over which this state will persist. By what mechanisms an open system that is clearly operating far from equilibrium will attain such an equilibrium state has never been defined. The wide range of significant time scales in Earth’s climate system, with significant interactions between sub-systems operating at different time scales, and the chaotic trajectories associated with all changes and interactions, present a challenging modeling and calculational problem. The dynamics make a big difference in the overall energy balance and are chaotic.
In computational fluid dynamics terms the problem with the boundary condition argument is that it ignores the influence of the dynamics of the system on the distribution of energy within the system. Turbulence for example is a powerful mechanism for energy redistribution from inertial scales to ultimately heat. The distribution of energy in the system can have a quite significant influence on the energy leaving or entering the system. For example turbulent processes such as convection and clouds play a large role in terms of the top of atmosphere energy balance. Another example is Rossby waves which modulate the equator to pole temperature gradient and thus the distribution of energy radiated to space. And this distribution determines the total amount radiated to space.
Computational fluid dynamics is based in large part on the insight that understanding the dynamics and thus the distribution of energy is critical to understanding the system even in total energy balance terms.
Recently, several papers have appeared that document some of the effects of dynamics on total energy balance. Zhao et al [Ming Zhao et al., 2016: Uncertainty in Model Climate Sensitivity Traced to Representations of Cumulus Precipitation Microphysics. J Climate, 29, 543-560] dealt with the effect of microphysics of cloud sub grid models on the ECS of a GCM and found that the effect was significant. Similarly, when the French IPSL modelling group recently improved the clouds and convective parameterization of its main model, the ECS reduced (per AR5 Table 9.5) from 4.1 C to 2.6 C. While the changes are small in terms of total energy fluxes, they are quite significant when compared to the total flux perturbations induced by changes in greenhouse gases and thus are quite important to what we care about in the system.
The dynamics make a big difference in the overall energy balance and are chaotic. Thus feedbacks mediated by these dynamics are critical to determining the effect of the forcings. Quantifying the dynamics is critical to the, often, small perturbations we are really interested in. Thus, the “boundary value” problem explanation is simply wrong and likely to sow confusion.
Here’s the Question
Within the framework of GCMs and all other factors summarized above,
- (1) identify the boundaries of the physical domain that are represented in GCMs, e.g define the surface B which constitutes the invariant boundary of the whole system under study.
- (2) describe the Boundary Condition(s) (BCs) that are imposed on the boundary, and
- (3) describe the equation(s) and dependent variable(s) to which each BC that applies at the boundary B, e.g prescribe explicitely N functions Gi(x,t) defined on B such as for all dynamical parameters Fi(x,t) defined on B we have Fi(x,t) = Gi(x,t) for all x included in B.
Suggest a mechanism allowing to obtain a unique solution of the dynamical equations depending only on the boundary conditions and justify the time scale T at which this mechanism would start to operate.
G. Browning and H.-O. Kreiss, (1985). Numerical problems connected with weather prediction. In: Murman E.M., Abarbanel S.S. (eds) Progress and Supercomputing in Computational Fluid Dynamics. Progress in Scientific Computing, Vol 6, p. 377-394. BirkhÃ¤user Boston
H.-O. Kreiss, (1980). Problems with different time scales for partial differential equations, Comm. Pure Appl. Math., Vol. 33, p. 399â€“439.
J. Oliger and A. SundstrÃ¶m, (1978). Theoretical and practical aspects of some initial-boundary value problems in fluid dynamics, SIAM J. Appl. Math., Vol. 35, p. 839â€“866.
R. Byron Bird, Warren E. Stewart, and Edwin N. Lightfoot, (1960). Transport Phenomena, John Wiley & Sons, inc. New York.
Marlow Anderson, (2009). Who Gave You the Epsilon? & Other Tales of Mathematical History, The Mathematical Association of America, Washington DC.