While mathematicians have often used bifurcation analyses and numerical continuation to study the behavior of nonlinear systems, the equations are often simple such as Duffing’s Equation or the Logistic Equations. Here we apply these techniques to an equation set very close to the single column models embedded in forecast and global climate models. This paper reports on an analysis of the coupling of the atmospheric boundary layer to the surface using techniques of nonlinear dynamical systems specifically numerical continuation. The results show that the equations support multiple solutions in certain parameter spaces and that solutions are highly sensitive to certain parameters that may not be specified with confidence. For the stable boundary layer when geostrophic speed is used as a bifurcation parameter, two stable equilibria are found--a warm solution corresponding to a high-wind regime where the surface layer of the atmosphere stays coupled to the outer layer, and a cold solution corresponding to the low-wind, decoupled case. Between the stable equilibria is an unstable region where multiple solutions exist. The bifurcation diagram is a classic S shape with the fold-back region showing the multiple solutions. This shape is a classical indication of a lack of predictability in a system. While steady state solutions for forecasts problems are often considered interesting but not necessarily relevant, our results show that time dependent solutions carried out within the parameter regions of multiple solutions show sensitive dependence on initial conditions and rather bizarre behavior. For example in the fold back parameter space time dependent solutions show temperatures may remain warm for most of the night then drop off in the early morning hours to the cold solution.

We have carried previous studies using a simple two-layer model of the atmosphere with a fairly complete surface energy budget. This allowed the dynamical analysis to be carried out on a coupled set of four ordinary differential equations. The present paper extends this work by examining additional bifurcation parameters and, more importantly, the analysis of a set of partial differential equations with full vertical dependence. Simple mathematical representations of classical problems in dynamical analysis often exhibit interesting behavior, such as multiple solutions, that is not retained in the behavior of more complete representations. In the present case the S-shaped bifurcation diagram remains with only slight variations from the two-layer model. For the parameter space in the fold-back region, the evolution of the boundary layer may be dramatically affected by the initial conditions at sunset. An eigenvalue analysis carried out to determine whether the system might support pure limit-cycle behavior showed that purely complex eigenvalues are not found. Thus, any cyclic behavior is likely to be transient.

In weather forecast and general circulation models the behavior of the atmospheric boundary layer especially the nocturnal boundary layer can be critically dependent on the magnitude of the effective model grid-scale surface heat capacity. Yet, this model parameter is uncertain both in its value and in its conceptual meaning for a model grid in heterogeneous conditions. Current methods for estimating the grid-scale heat capacity involve the areal/volume weighting of heat capacity (resistance) of various, often ill-defined, components. A bifurcation analysis with heat capacity as the bifurcation parameter shows that there are multiple solution regimes. Since heat capacity is a multiplying factor of the flux terms in the surface energy equations, it would appear that the model solutions would be highly dependent on this parameter. The bifurcation analysis shows that there are two regimes for sensitivity. For small values of heat capacity the system is highly sensitive to change in heat capacity. However, for large values of heat capacity the system shows little sensitivity to this parameter. This has important implications for numerical techniques for solving the surface energy budget. In most models only a very low order ODE solver is used for the surface temperature equation. This choice is probably rooted in the fact that since the flux terms are assumed to be uncertain a higher order scheme is not needed. However, the nonlinear analysis and numerical time dependent simulations show that in the small heat capacity regime numerical stability problems can exist.