# Combining Data and Models to Study Woody Plant Encroachment

Woody plant expansion into grasslands and savannas, which is accelerating worldwide, often affects ecosystem processes. Understanding and predicting the environmental and ecological impact of encroachment has led to a variety of methodologies for assessing its onset, transition, and stability. These methods generally rely on dynamical systems approaches.

We seek to understand the competition between alternate stable states of grasses, trees, or shrubs and the influence of climate, fire, precipitation, and livestock grazing. Of particular interest to us are climate and fire interactions occurring worldwide (both naturally and due to fire suppression), especially the impact of woody encroachment on prairie grasslands in the central U.S. The Konza Prairie Biological Station (KPBS) is a native tallgrass prairie preserve located in the Flint Hills of northeastern Kansas, a grassland region of steep slopes overlain with shallow limestone soil unsuitable for cultivation. The Flint Hills region encompasses over 1.6 million hectares and is the largest remaining area of unplowed tallgrass prairie in North America.

The KPBS is divided into various watersheds, each with distinct enforced conditions. These conditions include different rates at which each watershed is burned and possible grazing by cattle or bison. The station is an ideal place to blend data acquisition and utilization with models to understand phenomena such as woody encroachment.

Our research on woody encroachment combines both data and models, and our initial work builds upon the development of simple low-dimensional stochastically-forced models [4-5, 7-8]. A low-dimensional model with stochastic precipitation and fire disturbance can examine the complex interactions between precipitation and fire as mechanisms that may suppress or facilitate increases in woody cover [1].

To analyze the impact of fire and precipitation frequency and intensity, we employ local in-time Lyapunov exponents or so-called Steklov averages [2-3, 6] that assess convergence/divergence over different divergence timescales as a measure of relative stability or instability. Lyapunov exponents calculate the perturbation sensitivity of time-varying solutions of dynamical systems. By varying the parameters that control fire and precipitation, and comparing the dependence of quantities analogous to the largest Lyapunov exponent as a function of these parameters, we ascertain the relative control exerted on woody encroachment through these mechanisms.

Researchers have successfully employed the largest Lyapunov exponent in different areas of biology and ecology as an indicator of chaos (i.e., sensitivity to the initial state) and a measure of the relative sensitivity of parameterized systems. For example, Lyapunov exponents are indicative of a dynamical system’s exponential rate of change either away or toward a particular state of the system, such as a grass or woody state. A positive exponent illustrates an exponential rate away from the current state, i.e., an unstable state and transition to a stable one. A negative exponent indicates a trajectory toward the initial state, meaning that the state is stable. Therefore, if a given precipitation environment results in a negative Lyapunov exponent for woody fraction, the particular rainfall regime is stable for the expansion of woody species; a positive exponent signals an unstable regime.

Figure 1 displays the divergence time \(i\) versus the logarithm of divergence distance \(\tilde{y}(i)\)'s time-averaged logarithm, defined by

\[\tilde{y}(i) =< log \bigg(d_j(i)/d_j(0) \bigg)>\approx \lambda_1(i) \cdot i,\]

where the average is over \(j\) and \(d_j(i)\) denotes the divergence of the \(j\)th difference over a fixed length of time \(i\). \(\lambda_1(i)\) (analogous to the largest Lyapunov exponent) denotes the average divergence rate over time intervals of length \(i\). The algorithm in [6] allowed for computations using an appropriate delay coordinate embedding. Figure 1 illustrates the behavior of the normalized local Lyapunov exponents for a subset of the precipitation magnitudes, frequencies, and burn considered in [1] that corresponds to a 16-year fire frequency and high annual precipitation.

**Figure 1.**Normalized divergence as a function of the divergence time scale

*i*(in years) for parameter values that correspond to a 16-year fire frequency and high annual precipitation for grass assimilation, woody transpiration, top layer soil moisture, and bare soil evaporation. Figure credit: Nathaniel A. Brunsell and Erik S. Van Vleck.

Figure 1 also shows that these local type Lyapunov exponents produce an annual pattern with seasonal variations for some parameter configurations. Grass assimilation \((A)\), woody transpiration \((Tr_t)\), and top layer soil moisture \((\theta)\) have two peaks per year while other variables—such as bare soil evaporation \((E_s)\)—have one peak per year. Two peaks in these plots correspond to instability (sensitivity with respect to perturbations in that variable) over divergence times of zero-three months, six-nine months, 12-15 months, etc.; one peak per year corresponds to instability over divergence times of zero-six months, 12-18 months, etc. Perturbations are amplified by approximately \(e^{\tilde{y}(i)}\), so that they grow with time when the slope in Figure 1 is positive and decay when the slope is negative.

The results in [1] indicate that precipitation frequency is more significant than the intensity of individual precipitation events when controlling woody encroachment. Fire, however, has a much more dominant impact on limiting encroachment. These results suggest that fire management—in the form of more frequent burns—may be an effective strategy to slow the onset of woody species into grasslands. While climate change might predict a reduced potential for woody encroachment in the near future, these findings imply that a reduction in woody fraction may be unlikely when considering anthropogenic fire suppression.

Our next steps involve the incorporation of data into more complex land surface models (LSMs). Noah-MP is an LSM that employs multiple options for key land-atmosphere interaction processes. It contains a separate vegetation canopy defined by a canopy top and bottom; crown radius; and leaves with prescribed dimensions, orientation, density, and radiometric properties. The Noah-MP model prescribes both the horizontal and vertical density of vegetation using either ground- or satellite-based observations. It utilizes a dynamic vegetation model that allocates carbon to various parts of vegetation (leaf, stem, wood, and root) and soil carbon pools (fast and slow). The model distinguishes between C3 photosynthesis pathways (used by most plants) and C4 pathways (used more by semi-arid grasses), and defines vegetation-specific parameters for plant photosynthesis and respiration.

To understand transitions between these competing stable states, we are interested in quantifying several coupling metrics from Ameriflux data, remotely-sensed atmospheric profiles, atmospheric boundary layer height and soil moisture, and output from both a low-dimensional model and the Noah-MP model. These factors will help identify the coupling state and relationship to available energy and partitioning, vapor pressure, and temperature gradients to quantify each state’s stability and determine the underlying data and model sensitivity. Ultimately, this will increase our understanding of the physical mechanisms responsible for coupling dynamics in the central U.S.

**References**

[1] Brunsell, N.A., Van Vleck, E.S., Nosshi, M., Ratajczak, Z., & Nippert, J.B. (2017). Assessing the roles of fire frequency and precipitation in determining woody plant expansion in central U.S. grasslands. *J. Geophys. Res.: Biogeo., 122*, 2683-2698.

[2] Dieci, L., & Van Vleck, E.S. (2002). Lyapunov Spectral Intervals: Theory and Computation. *SIAM J. Numer. Anal., 40*, 516-542.

[3] Dieci, L., & Van Vleck, E.S. (2015). Lyapunov Exponents: Computation. In *Encyclopedia of Applied and Computational Mathematics* (pp. 834-838). Berlin Heidelberg: Springer-Verlag.

[4] Petrie, M.D., & Brunsell, N.A. (2011). The role of precipitation variability on the ecohydrology of grasslands. *Ecohydrol., 5*(3), 337-345.

[5] Porporato, A. (2003). Soil moisture and plant stress dynamics along the Kalahari precipitation gradient. *J. Geophys. Res., 108*(D3), 4127.

[6] Rosenstein, M.T., Collins, J.J., & De Luca, C.J. (1993). A practical method for calculating largest Lyapunov exponents from small data sets. *Phys. D: Nonlin. Phen., 65*(1-2), 117-134.

[7] Staver, A.C., Archibald, S., & Levin, S.A. (2011). The global extent and determinants of savanna and forest as alternative biome states. *Science, 334*(6053), 230-232.

[8] Xu, X., Medvigy, D., & Rodriguez-Iturbe, I. (2015). Relation between rainfall intensity and savanna tree abundance explained by water use strategies. *Proc. Nat. Acad. Sci., 112*(42), 12,992-996.