## Abstract

Stomatal resistance (*R*_{s}) forms a pivotal component of the surface energy budget and of the terrestrial biosphere–atmosphere interactions. Using a statistical–graphical technique, the *R*_{s}-related interactions between different atmospheric and physiological variables are resolved explicitly from observations made during the First ISLSCP (International Satellite Land Surface Climatology Project) Field Experiment (FIFE). A similar analysis was undertaken for the *R*_{s} parameterization schemes, as used in the present models. Three physiological schemes (the Ball–Woodrow–Berry, Kim and Verma, and Jacobs) and one operational Jarvis-type scheme were evaluated in terms of their ability to replicate the terrestrial biosphere–atmosphere interactions.

It was found that all of the *R*_{s} parameterization schemes have similar qualitative behavior for routine meteorological applications (without carbon assimilation). Compared to the observations, there was no significant difference found in employing either the relative humidity or the vapor pressure deficit as the humidity descriptor in the analysis. Overall, the relative humidity–based interactions were more linear than the vapor pressure deficit and hence could be considered more convenient in the scaling exercises. It was found that with high photosynthesis rates, all of the schemes had similar behavior. It was found with low assimilation rates, however, that the discrepancies and nonlinearity in the interactions, as well as the uncertainties, were exaggerated.

Introduction of CO_{2} into the analysis created a different dimension to the problem. It was found that for CO_{2}-based studies, the outcome had high uncertainty, as the interactions were nonlinear and the schemes could not converge onto a single interpretive scenario. This study highlights the secondary or indirect effects, and the interactions are crucial prior to evaluation of the climate and terrestrial biosphere–related changes even in the boundary layer perspective. Overall, it was found that direct and indirect effects could lead the system convergence toward different scenarios and have to be explicitly considered for atmospheric applications at all scales.

## Introduction

Experiments such as HAPEX-MOBILHY (Hydrological Atmospheric Pilot Experiment-Mobilisation du Bilan Hydrique) (André et al. 1986), FIFE (Sellers et al. 1988), EFEDA (ECHIVAL Field Experiment in a Desertification-threatened Area) (Bollé et al. 1993), and the Vegetation and Energy Balance Experiment (Niyogi et al. 1995; Raman et al. 1998) have asserted the significant impact that surface features, such as vegetation and soil moisture, have on the planetary boundary layer (PBL) structure. PBL processes are initiated (and modulated) through alterations in the surface thermohydraulic parameters. Hence, modeling efforts at all scales—micro (Su et al. 1996), meso (Alapaty et al. 1997), and global (Sellers et al. 1996; Randall et al. 1996)—attempt to realistically represent the responses of surface-induced changes (Jarvis 1976; Deardorff 1978; Avissar et al. 1985; Dickinson et al. 1986; Sellers et al. 1986, 1996; Wetzel and Chang 1988; Noilhan and Planton 1989; Xue et al. 1991; Acs 1994; Bosilovich and Sun 1995; Viterbo and Beljaars 1995; Alapaty et al. 1997; Niyogi et al. 1997a, among others). In addition to the meteorological applications for terrestrial biosphere–atmosphere interaction studies, physiologically intensive efforts are under way (Farquhar et al. 1980; Ball 1987; Finnigan and Raupach 1987; Meyers and Paw U 1987; Goudriaan 1988; Lynn and Carlson 1990;Grantz and Meinzer 1990; Raupach 1991; Collatz et al. 1991, 1992; Kim and Verma 1991; Baldocchi 1992, 1994; Dougherty et al. 1994; Jacobs 1994; Nikolov et al. 1995; Su et al. 1996; Sellers et al. 1996; Randall et al. 1996; Cox et al. 1996; Foley et al. 1996).

The meteorological and physiological treatments of the terrestrial biosphere–atmosphere interactions have certain subtle differences that need to be addressed. As shown in Fig. 1a, the generic pathway for the communication between the plant response and the environment is through the changes in the forcing that perturbs the basal state of the biotic system. The biota, in turn, attempt to dampen the impact the perturbation could have by adopting to the forcing themselves. Consequently, a combined (or effective) impact of the changes in the biotic signature and the external feedback is manifested as a plant response. The modulation attempted by the biota could provide a feed-forward signal to the forcing and thus modulate the environment. The effectiveness of the feed-forward signal would depend on various factors, such as the amplitude of the signal, the spatial coverage of the biota, and the modulation from other abiotic forcing, such as the soil moisture and radiation. In response to this communing, a rather straightforward coupling is postulated in the meteorological approach that is generalized in Fig. 1b. The example taken is apropos stomatal resistance (*R*_{s}). Here, *R*_{s} regulates the vegetation evapotranspiration, or the latent heat flux (LHF), thus regulating the surface energy budget and the overall PBL structure (Jacobs 1994; Alapaty et al. 1997; Niyogi 1996; Niyogi et al. 1996, 1997a). In the meteorological approach, referred to in Fig. 1b, the environmental response is accommodated in parameters such as temperature, radiation, and humidity (Jarvis 1976). These external forces brace with the plant’s internal or basal characteristic responses such as the minimum stomatal resistance (*R*_{s,min}) and LAI (leaf area index). The changes in the stomatal resistances are in turn reflected in the ambient temperature and humidity modulation due to the nexus between the biota and the surrounding air–boundary layer conductance. On the other hand, the physiological schemes estimate *R*_{s} through a noncausal photosynthesis–evapotranspiration balance (Cowan 1982; Farquhar and Sharkey 1982; Dougherty et al. 1994; Niyogi et al. 1997c). Hence, as shown in Fig. 1c, the environmental response is perceived by the “meteorological” parameters, such as temperature, humidity, and radiation (and carbon dioxide), and is then communicated to the physiological changes through photosynthesis (*A*_{n}) and related variables, such as leaf or intercellular carbon dioxide concentrations and respiration. These responses deal with the characteristic basal nature of the biota through responses, such as the LAI, and with the maximum potential for carbon assimilation and radiation interception (Farquhar and Sharkey 1982). Thus, in comparison with the meteorological schemes in which the ambient environmental characterization and the basal modulation result in the stomatal response in a coupled fashion, the physiological schemes postulate that the coupling is between the ambient environment and the plant’s instantaneous or dynamic response and the basal characteristics.

The manner in which *R*_{s} is parameterized (meteorological or physiological) could affect the boundary layer and climate predictions significantly. Recent findings from studies such as Sellers et al. (1996) confirm the possible variability in the general circulation model (GCM) predictions due to an introduction of a carbon assimilation (photosynthesis) pathway. Similarly, for PBL analysis, Niyogi and Raman (1997, henceforth Niyogi and Raman) compared predictions from four *R*_{s} schemes with observations during FIFE (Polley et al. 1992). Of the four schemes used, three were physiological [Ball–Woodrow–Berry (Ball 1987; Ball et al. 1987, henceforth Ball–Woodrow–Berry), Kim and Verma (1991, henceforth Kim and Verma), and Jacobs (1994, henceforth Jacobs)] and the fourth was meteorological [Jarvis-type approach, as in Noilhan and Planton (1989, henceforth Jarvis type)]. Predictions from the three physiological schemes were in accord with the observations, while the Jarvis-type scheme showed some discrepancy between predictions and observations, as well as a sluggish response to the atmospheric changes. Studies such as Su et al. (1996), Nikolov et al. (1995), Baldocchi and Rao (1995), Vogel et al. (1995), Jacobs (1994), and Collatz et al. (1991) also affirm that the physiological photosynthesis-based schemes simulate the observed stomatal behavior more closely. However, the physiological stomatal resistance schemes are “deceptively simple” (G. Farquhar 1996, personal communication) and require information regarding photosynthesis rates and cuticular carbon dioxide concentration, which is not routinely available for meteorological applications. Though such information could possibly be generated in the future using data for foliage nitrogen availability (D. Balddochi 1996, personal communication; Schulze et al. 1994; Niyogi et al. 1997), the procedure is still exploratory. Additionally, studies such as Jacobs or Niyogi and Raman show that the classical Jarvis-type approach for the meteorological scheme does not account for explicit interactions or secondary feedback processes in the system dynamics (compare Figs. 1b,c). Thus, the meteorological scheme could be argued to be deficient compared to the physiological approach. Our present knowledge and hence the interpretation of such interactions is limited and intuitive. There is an immediate need to study such atmospheric interactions apropos the two *R*_{s} approaches, as they affect interpretation of the model results and computational efficiency (Randall et al. 1996). This study attempts to resolve the interactions taking place within the physiological and meteorological–environmental parameters resulting in a stomatal response.

The following section describes the methodology applied to resolve interactions within the atmospheric system using a fractional factorial (FF) approach (Niyogi 1996; Niyogi et al. 1995, 1996, 1997a,b; Niyogi and Raman 1998, unpublished manuscripts). Section 3 discusses the interactions extracted using the FF design from the field observations made during FIFE (Polley et al. 1992). Section 4 describes the interactions mimicked by the four *Rs* parameterizations tested by Niyogi and Raman using the FIFE input parameters. Conclusions for this study are presented in section 5.

## Resolving interactions

Direct and secondary feedback or interactions always exist with a varying significance in the atmospheric system. Observations and measurement protocols mostly address the conformity to an effective response of both the direct and indirect effects. Thus, the direct effects are causatively perceivable under an assumption that indirect effects are less dominant. The possibility for this association stems from the ill-posed nature of the indirect effects or the secondary feedback–feed-forward coupling in the system. Understanding and representing the secondary feedbacks in one of the terrestrial biosphere–atmosphere links (stomatal response) is the focal point of this study.

The feedback loop can be simplified as a sum of 1) the main effect and 2) the indirect equilibrium sought with other system variables. Hence, if *P*_{1} is the variable to be studied in a system that has two other dominant variables—say, *P*_{2} and *P*_{3}—then the effective interactions can be expressed as

where the *k*_{1}(·) term is the main-effect term and the *k*_{2}(:) terms are the secondary feedback of the interaction (Box et al. 1978; Haaland 1989; Niyogi et al. 1996). It should be noted that although *k*_{1} and *k*_{2} are represented as two separate terms, their effects are intimately linked, latently coupled, and in principle cannot be isolated through a sensible measurement. Hence, a majority of the analyses performed combine the two terms together along with the direct noninteractive feedback.

Analysis of the dynamic system thus has been a reflection of the combination of the direct and indirect feedbacks, adopting a classic one-at-a-time (OAT) approach. [See Stein and Alpert (1993) or Box et al. (1978) for an elaboration.] In the OAT approach, *P*_{1} is varied and the response is recorded; then, resetting *P*_{1} to its original value, *P*_{2} is varied and the process is repeated. Upon viewing these responses, significance of the individual variables and causal relations between them are often implied. Haaland (1989) has demonstrated that by resolving interactions an enhanced understanding of the system “strategy” (feedback and feed-forward coupling) can be found. Similarly, within the atmospheric framework, recent studies (Stein and Alpert 1993; Alpert et al. 1995; Niyogi 1996; Niyogi et al. 1995, 1996, 1997a,b) have also shown the significance of resolving interactions. Alpert et al. (1995) conclude that without resolving interactions the results are often prone to be biased for the atmospheric system.

A crucial aspect in studying the soil vegetation–atmosphere transfer processes is understanding the stomatal resistance–related interactions. Stomatal response could represent a coupled interface for the biotic and abiotic components of the interactive transfer. This is because the stomatal response we perceive may be an equilibrium sought by the foliage between photosynthesis and transpiration as a strategic interactive balance with other stimuli, such as radiation or humidity (Farquhar and Sharkey 1982). Therefore, a change in any variable, say foliage temperature, may alter the equilibrium. It can be hypothesized that as a dynamic response to the change in the foliage temperature [quantified through main interaction or the *k*_{1} term in Eq. (1)], a secondary interactive change [inferred through the *k*_{2} term in Eq. (1)] will be manifested in another variable such as humidity, which itself has a *k*_{1}- and a *k*_{2}-type impact on the stomatal response. In an attempt to resolve and analyze such interactions for the stomatal behavior, data are considered from the observations of Polley et al. (1992) and from the model analysis of stomatal resistance for the Jarvis-type scheme and the three physiological schemes (Ball–Woodrow–Berry, Kim and Verma, and Jacobs) from Niyogi and Raman (Part I of this study).

The methodology followed is similar to the graphical–statistical approach developed by Niyogi et al. for the atmospheric framework (Niyogi 1996; Niyogi et al. 1995, 1996, 1997a,b,c; Niyogi and Raman 1998, unpublished manuscript). Using the FF approach, a high-resolution design (resolution V) is adopted (Box et al. 1978; Haaland 1989). In a resolution V design, all the main effects or direct interactions (*k*_{1} type) and the two-factor or indirect interactions (*P*_{1}:*P*_{2}, *P*_{2}:*P*_{3}, *k*_{2} type) are explicitly resolved along with statistically significant three-factor interactions (*P*_{1}:*P*_{2}:*P*_{3}, which are also *k*_{2} type), if present (Haaland 1989). Conforming to the design, combinations of the high and low variable settings are selected. For a combination, the resulting stomatal resistance is taken as the “effect” of the system. For this unsaturated design approach, the effect is subjected to main-effect analysis to quantify the impact due to *k*_{1}-type terms [see Eq. (1)] and to interaction explicit“Pareto” analysis (a representation with the variance ranked in a descending order) to resolve the *k*_{2}-type interaction terms (Haaland and O’Connel 1995; Niyogi 1996). The extracted interactions are further examined through interaction plots. The following section details the application of this approach for the observations made during FIFE.

## Interactions in observations

The observations used for the analysis were obtained as a part of FIFE (Sellers et al. 1988, 1992) and are documented by Polley et al. (1992). For different vegetation stages, observations of stomatal resistance and related factors, such as photosynthesis rate, CO_{2} concentrations, net radiation, temperature, and photosynthetically active radiation (PAR) were made over a C4 type of vegetation. Typically, gas exchange measurements over the foliage surface provided the physiological measurements, while the meteorological data were obtained from automated stations. [For details on the measurements, readers may refer to Polley et al. (1992) and Niyogi and Raman.] In this analysis, only the data for daytime observations were considered. To formulate a statistically robust hypothesis, diverse, yet representative, variable ranges are suggested (Box et al. 1978). To obtain the diverse combinations of variable settings from the limited observations, data from three different vegetation stages (days 157, 187, and 227) were compiled. Also, the observations were not separated for the three major generic families of the C4 vegetation type present [*Andropogon gerardii* Vitman, *Panicum virgatum* L., and *Sorghastrum nutans* (L.) Nash]. Studies made on the C4 vegetation in the FIFE study region (see Polley et al. 1992; Sellers et al. 1992, as well as the FIFE special issue in the *Journal of Geophysical Research,* 1992) show that the vegetation exhibits remarkably similar functional response to the atmospheric and biochemical changes, hence, such a grouping in a single dataset is acceptable. The data were checked for outliers and transformations, following Niyogi et al. (1997a), and then were used for developing combinations.

To resolve terrestrial biosphere–atmosphere interaction from observations, this study was divided into two subparts: meteorological and physiological. In the “meteorological” study, the variables considered were observations of stomatal resistance (*R*_{s}), ambient temperature (Temp), net radiation (Radn), and three humidity responses. The three humidity responses were vapor pressure deficit (VPD), relative humidity (RH), and ratio of vapor pressure deficit to maximum vapor pressure deficit (Drat). Here, the vapor pressure deficit was taken as the difference between the estimated saturated vapor pressure at the surface temperature and the vapor pressure at the ambient dry-bulb temperature, and it was used in the Jarvis-type scheme (Niyogi and Raman, Part I of this study). For the “physiological” study, the variables considered were observations of stomatal resistance (RsPhy), PAR, photosynthesis rate (Phot), leaf CO_{2} content (fCO_{2}), and foliage temperature (Tf). Relative humidity, estimated from the observations following a methodology given in Stull (1995), was considered as a humidity descriptor in the Ball–Woodrow–Berry model. Drat was taken as the humidity response in the Kim and Verma and also in the Jacobs models. [For maximum vapor pressure deficit, a value of 45 g kg^{−1} was considered, following Choudhury and Montieth (1986).] These variables are not routinely measured for meteorological studies and hence are considered physiological, as discussed in section 3b.

### Analysis of meteorological observations

The meteorological variables used for studying interactions were temperature (Temp, °C), relative humidity (RH, %), vapor pressure deficit (VPD, g kg^{−1}), radiation (Radn, W m^{−2}), and Drat. For these five variables, a resolution V design of the type FF0516 (fractional factorial design for five parameters and 16 combinations) as shown in Table 1 (Haaland 1989) was used. The “+” (high) and “−” (low) values for the variable setting required for generating the combinations consistent with the design could be obtained in a number of different ways. A simple way could be an intuitive or a subjective allocation of a value as being “high” or “low.” A more consistent approach, adopted in this study, is to segregate the data on the basis of a qualifier such as mean, median, or percentile (P. Haaland 1997, personal communication). The variable range in this study showed a lognormal distribution, and hence, the median was adopted as a delimiter for allocating high and low settings of the variable value. Thus, for example, the *R*_{s} value (in units of s m^{−1}) was taken as the effect for combination 1 (Table 1) and for the combination of low temperature, low relative humidity, low vapor pressure deficit, low radiation, and high Drat. On the other hand, the observed *R*_{s} for a condition corresponding to high temperature, high relative humidity, high vapor pressure deficit, low radiation, and low Drat was taken as the effect corresponding to combination 15. If more than one set of observations were found consistent with a combination, then a set for which the parameter values were farther away from the median was considered. The *R*_{s} observations corresponding to the design were then subjected to a graphical analysis through main effects and Pareto plots, as shown in Fig. 2.

In Fig. 2 the main-effects plot brings out the significance of the direct interactions [*k*_{1} term in Eq. (1)] of the parameter. The slope of the line joining the mean of the low and high variable setting gives the relative significance of the variable and the tendency of the system to respond to a change in the variable. As may be seen for radiation (Radn in Fig. 2), the slope of the line joining the effect corresponding to the lower and higher setting, respectively, indicates that radiation is a significant interaction for *R*_{s} as a response. Also, lower radiation generally results in a lower *R*_{s} main effect, whereas higher radiation increases the *R*_{s} main effect. In contrast to this direct relation, all of the other four variables (temperature, vapor pressure deficit, relative humidity, and Drat) are inversely related to the *R*_{s} main effect. All the three humidity descriptors have a comparable response, with RH being the strongest. The Pareto plot explicitly accounts for the secondary interactions in the system parameters. In the Pareto plot, terms such as Temp:Drat indicate interaction between Temp and Drat. For this case, the temperature and Drat interaction is the most sensitive or significant of the variables. The size of the effect quantified (as −58.2 for Temp:Drat or −57.4 for RH) from the analysis of variance approach (see Haaland 1989; Box et al. 1978, 314, for a discussion) attributes the importance to the parameter for the system effect (*R*_{s} in this study). Thus, in terms of the direct or main effect of the variables, the order of significance is RH, Radn, Drat, Temp, and VPD. Additionally, the interaction terms such as Temp:Drat, VPD:Radn, Temp:VPD, and Temp:RH are also quite significant. The two-factor interactions for the different parameters are shown in Figs. 3a–d.

Figure 3a shows the interactions for changes in the ambient temperature. Accordingly, Drat and Temp interact synergistically. That is, the impact of either Drat or Temp cannot be isolated for interpretation of the *R*_{s} outcome. Again, we see a similar response of the other two humidity descriptors: RH and VPD. At higher values of the parameter they become independent of temperature, while at lower ranges they are strongly temperature dependent in contributing toward the effective *R*_{s}. Conversely, for low temperature, both RH and VPD have a strong effect on *R*_{s}, which decreases with increasing temperature. On the other hand, Drat has a larger impact on *R*_{s} with increasing temperature. Figure 3b shows the interactions with Radn as the criterion. In the figure, as expected, VPD and Radn are strongly interactive. For lower vapor pressure deficit, *R*_{s} is strongly coupled to the changes in the Radn values. At higher VPD values, however, this VPD–Radn correlation is lowered. Thus, models with the VPD approach need to have considerably lower uncertainty in Radn parameterizations to reduce the errors associated with evapotranspiration estimates. Relative humidity and Temp show a similar behavior, with higher values being better correlated with Radn changes in governing *R*_{s} as a response. Drat appears to be a sensitive variable always showing Radn dependence. In Fig. 3c, the changes in *R*_{s} with changes in RH are delineated for the four other parameters. At lower RH, *R*_{s} is linked with Drat and Temp, while as RH increases, the VPD and Radn effects start dominating the *R*_{s} control. Once again, close interaction between Drat and Temp is observed. Also, except for VPD, at higher parameter settings, there is poor interaction with RH for effective *R*_{s}. That is, at higher RH values, the observed *R*_{s} is principally related to RH itself or VPD. The VPD dependence is shown in Fig. 3d. Drat, by definition, is directly related to VPD by a constant. Hence, no interaction is seen between the two variables. There are, however, strong interactions between Temp and VPD, and Radn and VPD. As also seen in Fig. 3c, at higher VPD there is relatively less impact of other variables on *R*_{s} for the RH analysis. Hence, both RH and VPD, consistent with their definition, are similar in their performance in the biospheric setup for generating *R*_{s} as a response. It thus appears from these results that the significance of humidity or vapor pressure deficit as a humidity descriptor is not largely different from the other (however, see also Aphalo and Jarvis 1991; Monteith 1995a,b).

This analysis for the meteorological variables provides some interesting insights into the temperature, radiation, and humidity descriptors with *R*_{s} as a response. The next subsection outlines the experiment with physiological observations.

### Analysis of physiological observations

For the physiological analysis, the variables considered are observations of PAR (*μ*einsteins m^{−2} s^{−1}), Phot (*μ*mol of CO_{2} m^{−2} s^{−1}), fCO_{2} (*μ*mol mol^{−1}), and Tf (°C). The response is *R*_{s} (the inverse of conductance that is in mol m^{−2} s^{−1}). (Stomatal resistance was analyzed in molar units so as to be consistent with the manner in which physiological data are analyzed.) For the four parameters, a resolution V design is of the form FF0416 (Haaland 1989). The design selected is shown in Table 2. As in the earlier case, the sample median is the criterion for qualifying a value as higher (+) or lower (−) from the observations available.

The main effect and Pareto plots for this analysis are shown in Fig. 4. As expected (Farquhar et al. 1980; Polley et al. 1992), photosynthesis rate and PAR have the maximum direct interactive response on *R*_{s}. Thus, from the interactions delineated in Fig. 4 with a lower PAR or photosynthesis rate, the resistance offered is higher, while with increasing foliage temperature and the intercellular (fCO_{2}) content, the observed resistances are higher. However, the Pareto plot indicates fCO_{2}, and, photosynthesis-rate interaction is significant. Hence, the fCO_{2} outcome cannot be interpreted as independent of assimilation, as it can in the main-effect analysis. In fact, significant fCO_{2} interactions are seen for all of the variables considered, with a significant Tf–fCO_{2}–Phot triple interaction term. This implies that fCO_{2} is an important variable for the stomatal resistance response. However, due to the interactive nature of fCO_{2}, model stability may be affected and a trade-off has to be achieved. From the Pareto analysis, both the increasing PAR and the photosynthesis rate are related to the decreased outcome of the *R*_{s} value. Higher Phot and fCO_{2} interaction as well as the fCO_{2} and Tf interaction are linked with reduced *R*_{s}, while the triple interaction (Phot–Tf–fCO_{2}) is positively correlated to *R*_{s}. This suggests that the majority of meteorological studies that attempt to present *R*_{s} analysis independent of fCO_{2} could have a strong bias due to the interactions observed.

Another interesting feature is the PAR and Phot interaction. Both these terms, independent in direct interaction are inversely related *R*_{s}, while the secondary interaction term is positively correlated. Hence, any model considering these two effects additively may have better performance if an interaction term of the functional form (PAR × Phot) is also included. This could be one of the perceivable reasons that a photosynthesis response term cannot be directly linked to the existing diagnostic Jarvis-type approach through a functional form, suggesting a need for an explicit photosynthesis equation in the model.

In the present analysis, with physiological variables, we could identify various significant indirect interactions. Figures 5a–d show the two-factor interaction plots obtained from the Pareto analysis (Fig. 4). In Fig. 5a the interactions of the physiological variables with PAR are delineated. From the main effect and the Pareto plot (Fig. 4), both PAR and fCO_{2} have an inverse relation with *R*_{s}. The two-factor interaction plot shows a synergistic interaction between PAR and fCO_{2}, indicating that their combined change will act more intensely than the variables will individually. Also, both Tf and PAR are closely linked and yet they do not interact directly. A similar independence is seen for PAR and *R*_{s} when the photosynthesis rate is higher. However, as PAR increases, the photosynthesis-rate interaction of *R*_{s} decreases relatively. The Tf–fCO_{2}–Phot three-way interaction was found to be significant for the observations considered. Also, the interaction was found to be antagonistic with the main-effect and the Pareto-plot interactions (Fig. 4), showing an opposite dependency. Figure 5b views this dependency in more detail. For low photosynthesis rate and low Tf, *R*_{s} is strongly linked with the fCO_{2} value. The fCO_{2} and Tf interaction is synergistic, and hence, the fCO_{2} dependence could also be important for higher Tf values. On the other hand, for a higher photosynthesis rate (Phot = +), Tf, and fCO_{2} do not interact as intensely, though they each directly influence the *R*_{s} values. A similar feature related to Tf and the inhibited photosynthesis rate has been reported by Jacobs (1994). Notice that the Tf and fCO_{2} interaction is antagonistic. Therefore, the inverse dependence (lower fCO_{2} or Tf leading to higher *R*_{s}) may be more intense for the lower-photosynthesis-rate case and may decrease considerably with an increasing photosynthesis rate. Figure 5c shows the mean Phot–fCO_{2} and Phot–Tf interaction, which suggests an intense coupling for *R*_{s} outcomes as a whole. For higher photosynthesis rates, it was observed that *R*_{s} tends to be independent of all other variables (fCO_{2} and Tf) except Phot, and the same is true for PAR. For lower Phot, *R*_{s} is linked with PAR, but this dependence decreases as the photosynthesis rate increases. Earlier studies (Farquhar et al. 1980; Collatz et al. 1991, for simulation outcomes; Polley et al. 1992, for observational analysis) suggest that *R*_{s} outcomes are strongly related to the photosynthesis rates. This conclusion might be appropriate only for the high photosynthesis rate. As with cases for lower photosynthesis rates, all of the physiological variables have a significant dependence on the *R*_{s} outcome. In other words, for regions with climatically low carbon assimilation (photosynthesis) rates, the indirect effects will assume higher significance. From the main-effect plot, fCO_{2} had an inverse direct relation with *R*_{s}. However, it is hypothesized that increasing CO_{2} may result in increasing *R*_{s} (Houghton et al. 1990; Jacobs 1994). This could be an interactive effect of the increased coupling of *R*_{s} with other variables, such as photosynthesis rate, Tf, and PAR. Thus, the impact of a CO_{2} rise might be related not only to the direct effects of CO_{2} changes but also to the secondary interactions coupled to the enhanced sensitivity with other parameters for higher CO_{2} availability. The overall response then may not be interpreted simply in terms of increasing the resistance but also in terms of an altered (higher or lower) stomatal response, depending on the secondary interactions. Hence, when dealing with simulations involving the doubling of CO_{2}, it might be necessary, in specifying the physiological parameters, to resort to higher degrees of sophistication other than the bulk approach as applied in the present GCM analyses (Paw U et al. 1996) or the terrestrial biosphere–atmosphere interaction parameterization.

This part of the analysis presents the interactions resolved for observed data during FIFE. The next aspect of this study describes the ability of the parameterization schemes to mimic the interactions extracted from the observations.

## Interactions mimed by different *R*_{s} schemes

For analyzing the interactions in the stomatal resistance parameterizations, the output from Part I (Niyogi and Raman) is used. Thus, analysis of the meteorological scheme (Jarvis type) and the three physiological schemes (Ball–Woodrow–Berry, Kim and Verma, and Jacobs) is done for resolving interactions.

### Meteorological schemes

In the Jarvis-type scheme the stomatal response is modeled principally as a function of the environmental parameters. The parameterization is of the form

where *R*_{s} is the stomatal resistance, *R*_{s,min} is the minimum stomatal resistance taken as a constant, Rn is the net radiation, *w*_{2} is the deep soil moisture availability, VPD is the vapor pressure deficit, and Temp is the ambient temperature. As compared to the other parameters, deep soil moisture can be assumed to remain constant for a time span of physiological interest (a few hours). Thus, along with *R*_{s,min}, *w*_{2} is treated as a constant in this analysis, and a design is generated for Rn, VPD, and Temp. The design selected is shown in Table 3 and is of the form FF0308, with resolution V (Haaland 1989). Using the Jarvis-type simulation results from Niyogi and Raman corresponding to each design input combination, an *R*_{s} value is selected as an effect. This resulting *R*_{s} effect is then used for the interaction analysis.

Figure 6 shows the main effect and the Pareto plot for *R*_{s} from the Jarvis-type-based outcome. Accordingly, for the data considered, Temp and VPD show a positive interaction with *R*_{s}, while Radn is negatively related to *R*_{s}. In Part I of this study, we identified that the Jarvis-type scheme could not simulate the observations over the C4 type of vegetation under low-photosynthesis-rate scenarios. Therefore, it is important to examine the interaction analysis results. Interestingly, compared to the observations (Fig. 2), none of the parameters agrees in terms of the interactions, as the parameters show opposing behavior. From the Pareto plot for the *R*_{s} outcomes, it can be seen that the Jarvis-type scheme treats VPD and Radn as the crucial parameters. This is in agreement with the earlier results of Niyogi et al. (1996) in which the radiation-related terms were significantly interactive for the Jarvis-type scheme. It is also interesting to note that all of the main-effect terms show an opposite behavior for the model tendency as compared to the observations but that the interaction terms related to VPD show like behavior. Notice for the observations, however, that the interactions are antagonistic with the main effect and that the secondary interactions show the opposite tendency. For the Jarvis-type scheme, the interactions are synergistic. Niyogi and Raman found that the performance of the Jarvis-type scheme was poor in terms of simulating midday stomatal resistance over a C4 vegetation (although the performance was fair in predicting the ensemble mean). As the present analysis reveals, the incongruence in the tendencies for the scheme could be one of the causes of the poor feedback (Jacobs 1994).

In Figs. 7a–c, the two-factor interactions mimicked by the Jarvis-type scheme are obtained. Figure 7a shows the Temp–*R*_{s} interactions with VPD and Radn. For increasing Temp, the importance of VPD and Radn increases in this scheme. Comparison to the interaction plots obtained for the observations (Fig. 3) for high Temp cases shows that the model performance is different from the observations. The Temp–VPD–Radn interaction for the *R*_{s} outcome is further analyzed in Fig. 7b. For lower Radn, VPD affects *R*_{s} independent of Temp. The Temp dependence increases as VPD increases, while for increasing Radn, higher VPD is directly related to the *R*_{s} changes independent of Temp. Figure 7c details the VPD dependence on Radn and Temp. For the model outcome, VPD and Radn do not seem to interact, although they have a direct impact on the predicted *R*_{s}. The observations, however, show somewhat strong interactions between VPD and Radn. This interaction decreases with increasing VPD. Also, for the observations at higher Temp, *R*_{s} is not directly dependent on VPD, as is seen for the Jarvis-type outcome.

Thus, the Jarvis-type outcomes suggest that overall there is poor agreement with the observations in terms of the *interactive* response. It is expected that by understanding the behavior of the physiological schemes one can interpret a physiological-interactive response in the Jarvis-type analysis as an error or uncertainty term. However, our present understanding of the response of the parameters within the physiological schemes is not sufficient for providing information on interactions implicitly parameterized for a diagnostic prescription of *R*_{s}. The following section attempts to resolve the interactions for the physiological schemes. The outputs obtained from different models (Niyogi and Raman) are treated equivalently to observations in the earlier part of this study. The schemes examined are Ball–Woodrow–Berry, Kim and Verma, and Jacobs. The three schemes are assumed to surrogate the physiological approach of the stomatal response parameterization.

### Physiological approach

The physiological schemes for the *R*_{s} parameterization are governed by the photosynthetic processes. The physiological approach paired with the appropriate photosynthesis parameterization (Farquhar et al. 1980; Collatz et al. 1992) is known to simulate *R*_{s} outcomes close to observations (Dougherty et al. 1994; Vogel et al. 1996). Even in Niyogi and Raman, the physiological schemes performed well and followed observations closely. In this section, we will analyze the interactions as mimicked by the physiological schemes. For all three schemes—Ball–Woodrow–Berry, Kim and Verma, and Jacobs—the following four parameters are considered: a humidity descriptor, foliage temperature, net radiation, and foliage CO_{2} content. The design adopted for this analysis is of the form mf0412 (mixed fractional factorial analysis for four parameters with 12 runs). The mixed factor was considered since the humidity response for the parameter showed a distinct nonlinear distribution that could not be treated in a linear two-level system (unless normalized). The advantage of this design is that for one parameter three parameter levels can be selected (“high” and “low,” as in FF designs, with an additional “medium” parameter setting). The design used, following Haaland (1989), is shown in Table 4. The high (+) and low (−) values are again based on the sample median, as described for the previous analysis (section 3), and the intermediate (0) values are taken from the parameter values that are close to the sample median.

We will view the interactions for the various parameters in the three schemes, starting with the Ball–Woodrow–Berry scheme, followed by the Kim and Verma scheme, and finishing with the Jacobs scheme.

#### Ball–Woodrow–Berry

The humidity descriptor for this scheme is relative humidity (RH). The final equation for calculating *R*_{s} (inverse of *g*_{s}) is as follows:

where *g*_{s} is the stomatal conductance (the inverse of which is taken as *R*_{s}), *A*_{n} is the net CO_{2} assimilation or the photosynthesis rate, rh is the relative humidity at the leaf surface, and *C*_{s} is the CO_{2} concentration at the leaf surface. The remaining terms *m* and *g*_{o} are constants based on linear gas-exchange considerations (Ball 1987). The scalar concentrations at the leaf surface are often established as described in Nikolov et al. (1995) and Su et al. (1996).

The parameters selected for this analysis were RH, Radn, Temp, and CO_{2}. Corresponding to the input parameter combinations, as required for the mf0412 design (mixed factorial design with four parameters and 12 combinations) (Table 4), the *R*_{s} values predicted by the Ball–Woodrow–Berry scheme were selected. The main effect and the Pareto analysis outcome are shown in Fig. 8a. Figures 8b and 8c show the responses for the Kim and Verma model and the Jacobs model, which are discussed in the following sections. From Fig. 8a (the main-effect analysis) it is shown that RH and Radn are inversely related, while Temp and CO_{2} are directly related to the predicted *R*_{s}. The observed Temp and Radn effects are opposite to those seen for the Ball–Woodrow–Berry model. This scheme, however, had the best relative performance of the four schemes considered by Niyogi and Raman. In this analysis, though, only the RH direct interactions match the observations. From the Pareto analysis shown in Fig. 8a, the accurate estimation of foliage CO_{2} concentration is pivotal for this model scheme. In the observations, the direct effect of CO_{2} was to reduce *R*_{s}, while through indirect interactions the net effect was to increase the resistance effectively. The Ball–Woodrow–Berry scheme adopts this feature in a gross sense by directly relating the CO_{2} rise with increased resistance. In the Pareto plots, the terms with the L extension, such as RH.L or RH.L CO_{2}, indicate linear dependency, while terms with the Q extension, such as RH.Q and RH.Q Radn, indicate quadratic relations for the indirect interactions. Thus, the linear secondary interaction between RH and CO_{2} is prominent alongside the RH indirect effects. The linear and quadratic RH main-effect terms are inversely related to *R*_{s}. Both the L and Q RH–Radn interaction terms are directly linked with *R*_{s}. The RH–Radn interactions, however, are antagonistic, and reduced sensitivity of these terms in the scheme can make the parameterization robust. The scheme shows a different interaction response for linear and quadratic effects of Temp and RH. An interaction of the RH–Temp form for this scheme causes an increased *R*_{s}, while a higher-order interaction has the effect of reducing the resistance. This indicates for a nonlinear distribution of RH change that *R*_{s} could be overpredicted under low RH conditions and underpredicted under high RH conditions. This is an interesting feature. Evapotranspiration calculations consider the total resistance (the sum of *R*_{s} and the aerodynamic resistance *R*_{a}). Studies such as Collatz et al. (1991) suggest that under high RH environments the Ball–Woodrow–Berry approach could induce excess *R*_{a} (considering *R*_{a} relates with the boundary layer conductance). Hence, results from studies such as Collatz et al. (1991) suggest a possible antagonistic interaction between *R*_{a} and *R*_{s} for open C4 pathways in afternoon conditions. This implies that for evapotranspiration estimation the *R*_{s} underprediction and the *R*_{a} overprediction could reduce the error in the net sense. This feature needs to be verified from additional observations.

Figures 9a and 9b show the interaction plots related to the Ball–Woodrow–Berry outcomes. Consider the Temp–*R*_{s} interactions for this scheme. For higher Radn, the effect of change in Temp is insignificant, while with decreasing Radn, the temperature changes start affecting the stomatal response. This implies that the system attains stationarity under high radiation (midafternoon) availability. This stationarity is not disturbed with temperature changes. However, under low-radiation conditions (morning or early evening), temperature changes, such as from cumulus activity or advection, can be important sources of error in the evapotranspirative flux estimation in the boundary layer modeling. It follows that for higher Temp, Radn is still a significant parameter. In the Pareto analysis, RH was linked inversely with *R*_{s} as both a linear and a quadratic effect. We will analyze the RH interactions with three other model parameters considered in the design. According to the Ball–Woodrow–Berry parameterization, the sensitivity of the CO_{2} change is larger for lower RH. This sensitivity reduces with increasing RH. A similar interaction is noted for RH and Radn. For lower Radn and lower RH, resulting *R*_{s} values are high, and they decrease with increasing Radn *or* RH. Interestingly though, for increasing Radn *and* increasing RH, there is an initial decrease followed by an increased *R*_{s}. Overall, for the CO_{2} interactions, the Radn impact is lower for higher RH. For the Radn–Temp interaction we found that at lower Radn the Temp changes affect the *R*_{s} outcome. There is a strong interactive coupling in the RH and Temp. For low Temp, an increase in the RH will increase *R*_{s} until a certain threshold and then start reducing the resistance again. However, for high Temp, increasing RH persistently decreases the predicted *R*_{s}. This occurs until a limiting high RH value is reached, beyond which the Temp changes do not affect the *R*_{s} outcome.

Thus, for the Ball–Woodrow–Berry scheme, RH dominates the *R*_{s} prediction. Accordingly, the significance associated with a change in other parameters is inversely related to the RH availability of the vegetation environment. However, the results need to be qualified, as the indirect results are not well documented in prior studies and more direct observations need to be taken.

#### Kim and Verma and Jacobs

The other two schemes we will analyze for interactions are those based on the work by Kim and Verma, and Jacobs. In the Kim and Verma scheme, the *R*_{s} equation has the following form:

while in the Jacobs scheme the equation for *R*_{s} is

In Eqs. (3) and (4) *A*_{n} is the net CO_{2} assimilation and *C*_{s} and *C*_{i} are the CO_{2} concentrations at the leaf surface and within the leaf, respectively.

The parameters chosen for both schemes were Drat (*D/D*_{max}), Radn, and CO_{2}. Once again, the design used is the same as that for the Ball–Woodrow–Berry analysis shown in Table 4 (mf0412 type). The predicted *R*_{s} corresponding to the input combination and to the design is taken from Niyogi and Raman. Although the predicted *R*_{s} from the two schemes are different (generally the Jacobs scheme had explicit interactions with higher *R*_{s} outcome and larger scatter), they show a close similarity in the response to a parameter change in the main effect, the Pareto analysis, and the interaction plots. This similarity is not surprising since the two schemes stem from a similar approach based on the Drat control of *R*_{s}. Hence, in this section the two schemes are combined as Drat schemes rather than as a Kim and Verma scheme or a Jacobs scheme.

Figures 8b,c show the main effect and the Pareto plots for the two Drat schemes. For the Ball–Woodrow–Berry scheme, the vapor descriptor, RH, had a persistent inverse relationship with *R*_{s}. The Drat–*R*_{s} relationship is more nonlinear. Until a certain threshold intermediate value of Drat is reached, *R*_{s} increases and then further rises in Drat, resulting in a decreasing *R*_{s}. Temp is directly related to *R*_{s}, and so are Radn and CO_{2}. All of these features agree with the observations analyzed earlier.

From the Pareto plot, the quadratic Drat term and CO_{2} gather maximum significance. Interestingly the quadratic Drat term (Drat.Q) is negatively related to *R*_{s}, while the linear Drat term is positively related to *R*_{s}. A similar feature was seen for the Ball–Woodrow–Berry scheme, with its humidity parameter (RH) in the nexus with *R*_{s}. Drat is interactive with other parameters, as in the Ball–Woodrow–Berry scheme, where RH showed interactions with every other parameter. Considering the Kim and Verma scheme and the Jacobs scheme individually, the only difference in terms of the interaction response generated by the two schemes is through the Drat.Q–Temp term. In the Kim and Verma scheme, the quadratic Drat and Temp interaction increases *R*_{s}, while for the Jacobs scheme the interaction reduces the *R*_{s} outcome.

Figures 10a,b and 11a,b summarize the interactions for the Drat schemes of Kim and Verma, and Jacobs. Considering the Radn–Temp interaction, we found that as a main effect (Figs. 8b,c), both Radn and Temp directly increase *R*_{s}, while the interaction in the Pareto analysis revealed a negative tendency. The two-factor interaction plot shows that the tendency for the two terms is to generate an antagonistic and a proportional response. Thus, the net impact of Radn and Temp change may not be as pronounced as perceived by earlier studies (Kim and Verma, and Jacobs). Some differences are seen for the Ball–Woodrow–Berry scheme and the Drat approach. For the Ball–Woodrow–Berry scheme we found that for high Radn *R*_{s} is insensitive to Temp changes. In the Drat scheme, however, Temp is *always* a sensitive parameter. The Drat changes are nonlinear as compared to the linear RH dependence in the Ball–Woodrow–Berry scheme. In the Ball–Woodrow–Berry scheme, a strong RH–Temp coupling was noticed that had a combined impact on the *R*_{s} outcome. In the Drat schemes, the Drat–Temp interaction is limited, with low temperatures being associated with low *R*_{s}. With a limiting rise in Drat, *R*_{s} increases, and further rise in Drat reduces the resistance. The Temp sensitivity, however, persistently increases with increasing Drat, while the Radn–Drat interaction is more interactive. For lower Radn and increasing Drat, the predicted *R*_{s} also increases. The *R*_{s} increase is more rapid from intermediate to higher Drat values. For higher Radn with intermediate Drat values, *R*_{s} increases but reduces steeply for higher Drat. Overall, the Drat sensitivity is higher for higher Radn values. For the CO_{2}–Drat coupling, a higher CO_{2} results in an increased *R*_{s} value. Also with increased Drat the CO_{2} sensitivity increases. Both these aspects are consistent for the physiological schemes. However, CO_{2}–Drat coupling is essentially interactive and must be considered in the Drat schemes. For low CO_{2} and a change in Drat from low to intermediate values, *R*_{s} increases. With a further rise in the Drat values, the *R*_{s} values decrease rapidly. For high CO_{2}, a similar high *R*_{s} response is seen with the intermediate Drat values. However, further rise in Drat does not alter the *R*_{s} outcome significantly.

## Conclusions

Stomatal resistance (*R*_{s}) is one of the most significant terms in a PBL analysis over a vegetated surface. The coupled role of terrestrial biospheric–atmospheric interactions for *R*_{s} estimation was examined. Using a fractional-factorial (FF) approach, *R*_{s} observations during FIFE were used for evaluating the interactions in the meteorological and physiological parameters. A similar analysis was undertaken for different *R*_{s} parameterizations, such as the Jarvis-type meteorological approach and the RH and VPD humidity descriptor approach, as in the Ball–Woodrow–Berry model and the Kim and Verma (1991) and Jacobs (1994) models. Several significant interactions were resolved, and the implications of these findings are discussed in this section.

For meteorological applications (PBL studies, mesoscale circulation studies, weather forecasting) all of the *R*_{s} estimation approaches have qualitatively a similar interactive behavior. Thus, there is no significant, apparent difference in employing Drat, VPD, or RH as a humidity descriptor in the analysis. The Drat schemes are strongly coupled with temperature changes. This Drat–Temp interaction can make *R*_{s} predictions from such schemes prone to larger scatter and thus to higher uncertainty during scaling from leaf to canopy to landscape. For such scaling, RH (as employed in the Ball–Woodrow–Berry scheme) appears to be well suited, as nonlinearities are eliminated in the interactions. Because of this, the Ball–Woodrow–Berry scheme may be better posed for meteorological modeling. However, if the Ball–Woodrow–Berry scheme is applied in meteorological studies, it will be necessary to have a noniterative photosynthesis coupling. Such a diagnostic coupling will require explicit consideration for higher-order interactions between RH and radiation, and it may have its best performance under high assimilation rates. It is, however, interesting to note that the performance of the Jarvis-type scheme is also reasonably good for high photosynthesis conditions (Niyogi and Raman 1997). Under low assimilation rates various parameters were seen to interact simultaneously. This is an important aspect that needs to be studied further. It is possible that low assimilation rate is typical for some geographical domains such as the Tropics (Schulze et al. 1994). For such regions or periods, mesoscale assimilation schemes (Mahfouf 1991; Niyogi et al. 1996c) could propagate large uncertainties embedded in the model predictions. GCM results have to be interpreted cautiously for such domains if secondary interactions are not explicitly resolved.

In contrast to the meteorological applications, the implications are different for studies that consider the impact of CO_{2} changes. The interactions examined suggest that the model predictions of the altered CO_{2} scenario will depend on the approach one adopts (viz., RH or Drat based in the ecological models, or VPD as in the GCMs). In fact, carbon dioxide–related terms are significantly interactive and appear difficult to simulate or parameterize since such effects would tend to be chaotic and would not necessarily replicate themselves (Niyogi and Raman 1998, unpublished manuscript). The interaction effect can be partly simulated through an iterative approach, but this will further affect the model stability for diverse scenarios. This suggests that an analytical or a diagnostic approach, as required for weather forecasting and climate models, could have large inherent uncertainties with CO_{2} changes. Thus, though much progress has been made related to terrestrial biosphere–atmosphere interaction research, questions remain regarding the coupled and interactive effects of CO_{2} changes, which this study has identified as a dominant forcing in the response in the terrestrial biosphere.

## Acknowledgments

We would like to acknowledge stimulating interactions with Dr. Dennis Baldocchi of ATDD, ORNL, and NOAA; Dr. Graham Farquhar of the Australian National University; Dr. Peter Harley of NCAR; Dr. Russ Monson at the University of Colorado at Boulder; Dr. Jim Randerson at the Carnegie Institute of Washington; Dr. Hong Bing Su from the University of California, Davis; and Dr. A. Prabhu at the Indian Institute of Science. In particular, we would like to thank Dennis Baldocchi and A. Prabhu for reviewing the manuscript and for their detailed comments. Thanks are due to Kelly Funk for her editorial assistance and to Dr. Perry Haaland for his suggestions regarding the fractional factorial technique. We thank Dr. John Norman, Univeristy of Wisconsin—Madison, for his e-mail correspondence regarding the FIFE observations. The first author (DSN) would like to acknowledge the opportunity, provided by the University Center for Atmospheric Research, to participate in the Advanced Study Program on Terrestrial Ecosystem and Atmosphere (TEA) at the National Center for Atmospheric Research (funded by the National Science Foundation). The interdisciplinary discussions concerning physiological interactions and representations with the TEA participants provided useful insight into the physiological schemes and the atmospheric model coupling mechanism, and is warmly acknowledged. The FIFE data used in this study were obtained from the Biogeochemical Information Ordering Management Environment (BIOME) at the ORNL’s Distributed Active Archive Center (DAAC). This research was funded by the Division of Atmospheric Sciences of the National Science Foundation, the Naval Research Laboratory, and the North Carolina Super Computing Center.

## REFERENCES

_{2}, and sensible heat flux densities and scalar profiles over and within a soybean canopy. Bound.-Layer Meteor., 61, 113–144.

_{2}assimilation in leaves of C3 species. Planta, 149, 78–90.

_{2}enrichment on regional transpiration. Ph.D. thesis, Wageningen Agricultural University, Wageningen, the Netherlands, 179 pp. [Available from Wageningen Agricultural University, Duivendaal 2, 6701 AP Wageningen, The Netherlands.]

_{3}plants. Ecol. Modeling, 80, 205–235.

_{2}concentration. Preprints. 12th Conf. on Biometeorology and Aerobiology, Atlanta, GA, Amer. Meteor. Soc., 58–59.

## Footnotes

*Corresponding author address:* Dr. Sethu Raman, Dept. of Marine, Earth, and Atmospheric Sciences, North Carolina State University, Box 8208, Raleigh, NC 27695-8208.