Using models
to help design the optimal strategy of deficit irrigation on nut crops
Vincent P. Gutschick
Global Change Consulting Consortium
4904 Calabazilla Rd.
Las Cruces, NM 88011
Phone 5755712269
Email vince@gcconsortium.com
Website http://gcconsortium.com
Dates of
preparation: Jan./Feb. 2009
Prepared on
consulting contract to Prof. T. W. "Ted" Sammis, Plant and
Environmental Sciences,
Table of
contents:
A prologue on what we expect, and how and why we are combining
models and experiments
The role of models, or, of good models
Making models
robust and mechanistic
Model role 1:
interpreting experiments
Model role 2:
designing the experiments
Finding
the most fundamental descriptors of plant responses
Complex
models vs. simple models
Example of
final use  a first estimate of optimal deficit irrigation
Here
are the premises of the simple optimization model
From
these premises, we can compute…
Finally,
we can find the optimal strategy of deficit irrigation
We will now
make an estimate, putting some quantitative expressions into the premises
Details of
Objective (1): Develop the models
Model canopy
photosynthetic rates (A), transpiration rate (T), soil evaporation rate
(E), and , canopy N content
The
need for both complex and simple models
The importance of canopy structure.
Photosynthesis
by the whole canopy
The
complete environment of the leaf determines photosynthesis
Solving all
these coupled equations
The
canopy affects its own environment
Some
processes are still left out
Programming,
and accounting for nutfill
Using
only one time of day as representative
The
model allows the user to enforce any desired cutback in water use
Using
the models, "forward" and "inverse"modes
Inputs to
the complex model: sources of data, quality control, formatting, and units
A
note on potential changes in required inputs
Meteorological
and other environmental information
Physiology
and other leaf traits
Controls
over desired accuracy and over reuse of information
Inputs to the simpler, bigleaf model
Basic
structure of the code for
the complex model
Interpreting
and analyzing the output
The
quantities provided in the output of the complex model, and their units
The quantities provided in the output of the simpler model,
and their units
Details of
Objective (2): Verify models with new and past field measurements / design
the best experimental protocols
Getting best estimates of both the key input
parameters and the detailed stress responses
Important
soil variables to relate to stress and physiological responses to stress
Some other
parameters, less variable
Observing
dynamic stress responses
Scaling up
the measurements of stress responses to whole canopies.
A brief sketch of objective
(3): Use the models to parametrize an initial optimization model
How the
models need to be modified  a quick outline
Training in
the use of the models and use of equipment
The research project sketched below will rely on a synergy between modelling and field measurements. This document is intended to describe the modelling effort, which is being used to synthesize our understanding, to help design experiments, and to contribute heavily to making a final decisionsupport system. It includes an overview of the specific project, a primer on the use of modelling (less familiar to most of the team members than is field experimentation), and some important features of the models that I am developing for the project.
Nut
production from pecans, almonds, and pistachios figures heavily in the
economies of
Irrigation at less than the levels that mazimize yield and/or profits will have to be done. It remains possible to set the timing and amounts of irrigation in such a way as to maximize the yield within the seasontotal constraints on water use. At the same time, the future yield capacity can be preserved and the death of trees prevented. To develop these optimal schedules for a given fractional availability of water, we must understand how the trees respond to deficit irrigation and its detailed scheduling. We must get quantitative estimates of how water stress changes tree photosynthesis, its partitioning (to nutfill, maintenance respiration, net growth, and reserves), and the dynamics of N in leaves, soil, and reserves.
Consequently, 13 investigators in the 3 states obtained federal funding to use their research skills for all aspects of this project. The work includes field measurements of physiological responses to water and N treatments, development of monitoring methods, and modelling so that the limited experimental results can be synthesized into an optimization plan applicable in varying climates and soils, across interannual variations, and with varied amounts of water cutbacks.
Ted Sammis engaged me, through my consulting firm (see http://gcconsortium.com) to do the modelling, and to apply it further in helping to develop the best experimental design and to train instrumentation users so that they might obtain the most useful and accurate data.
Develop optimal schedules (timing and amounts) of irrigation and N fertilization that maximize
yield when irrigation water is cut to 50% (or other specified fraction) of normal.
This goal can be achieved if we formulate a whole set of objectives in modelling and field measurements. These tasks will be presented later. First, we need to develop some context  what are the important tree responses, what can models do to help us put together the best strategy from our complete picture of tree responses:
A
prologue on what we expect, and
how and why we are combining models and experiments.
Expected tree responses. If trees or other plants are given reduced water supplies, many physiological acclimations can come into play. One that starts quickly and persists is a reduction in stomatal conductance, g_{s}. This cuts leaf transpiration almost in proportion  not quite as much, because leaf cooling is reduced, and the rise in temperature raises the leaftoair gradient in watervapor pressure. The reduction in g_{s} also cuts leaf photosynthesis, but considerably less than proportionally  the stomatal resistance (inverse of conductance) is a much smaller part of the total pathway resistance for incoming CO_{2}. Consequently, wateruse efficiency (WUE), as the ratio of photosynthetic rate to transpiration rate, rises. There are a number of ways to view this, conceptually, mathematically, and in measurement. These aren't detailed here. This gives us the hope that a 50% reduction in water use need not lead to as large a drop in yield.
Another acclimation is a set of changes in allocation of growth (allocation of photosynthesized sugars, or photosynthate, and all the metabolic products dervied from these). Plants may divert more photosynthate to reproduction. Perennials may devote less to reserves, which has consequences for the next season's growth and reproductive yield. The allocation of growth and also of nitrogen to leaves and within their photosynthetic machinery may change.
The role of models, or, of good models. Models joined with experiments are an excellent way to synthesize what we learn in experiments and then to estimate the best irrigation strategies. Experiments alone are insufficient and inefficient. Experiments to induce tree responses to water stress are difficult, expensive, and risky  using many replicate trees means using a large area, and it entails a risk of longlasting damage or death . Furthermore, we need to cover a wide range of climates, interannual variations in weather, soil types, etc. A multifactorial experiment would be wholly impossible.
Making models robust and mechanistic.
Models that incorporate robust patterns of plant responses to the environment
let us escape this conundrum, with much lowered risk of misleading
ourselves. This statement implies that
the models must be mechanistic, not statistical. By robust, I mean response
patterns that have been shown to be common among different species and
conditions. One very strong example is
the relation of leaf photosynthetic rate to CO_{2} concentration
(partial pressure) inside the leaf and the kinetics of Rubisco enzyme (or, in
lower light, a series of photochemical steps all coming down to one parameter,
an electron transport capacity). The
famed Farquhar  von Caemmerer 
The key is to use the robust patterns, not patterns that appear only under a limited range of conditions. Once we have a model using the most robust set of process descriptions, we can apply the model in diverse conditions, often over a far greater range than those in which we performed our specific experiments (but we should continue to test the model!).
Model role 1: interpreting experiments. Models consequently work with experiments in two ways. First, they interpret the experiments, to extract patterns in the large. In our case, we hope that they tell us how nut trees respond to an arbitrary sequence of daily environments and an arbitrary sequence of management decisions, esp. on irrigation. They then can be used to find the optimal management responses, using the mathematics of optimization (a rich and rewarding subject). They can also be used to make guidelines for realtime responses by growers, that is, decisionsupport systems.
Model role 2: designing the experiments. We want the experiments to do two things  tell us if our descriptions of plant processes are correct, and tell us the numerical values of the physiological parameters that we will use in models to estimate responses in arbitrary conditions. In the former part, we are basically testing hypotheses. We can only do a limited number of tests, and here is where models show another great benefit. They help us design the experiments so that they let us test the most critical plant responses. Provided that our model is reasonably good (inclusive of important processes and accurate enough in their description), it reveals the processes  and the parameters that describe these processes  that have the greatest effect on our final outcome. Models can be run with changes in the process descriptions and we can see the predicted sensitivity of, say, nut yield to these changes. For example, these simulations may reveal that stomatal control parameters are less important to final nut yield than is leaf photosynthetic capacity, hence leaf N level. They could further reveal that some parameters, such as tree spacing have very little effect. (This is just a notional example. There are reasons to suspect tree spacing has a significant influence, because canopies must be closed for best growth, and closing canopies with large spacing means letting trees grow large individually. In this condition, they put a greater fraction of growth into supporting trunk and branches and less into nuts). To conclude on this point, we find the top few parameters and choose an experimental design that will best highlight their effects.
Finding the most fundamental descriptors of plant responses. In the second part, we are trying to identify the most fundamental, unchanging aspects of plant performance. A good model makes us think hard about what processes occur and what are the most basic descriptions of them. Consider transpiration. What is the best way to describe it, in terms of basic plant properties? Certainly, stomatal conductance is a most direct controller of transpiration, but g_{s} varies widely with the diurnal variations of the environment and on different leaves. We can't use it as a basic parameter. Consider, then, the BallBerry equation. If this truly captures most of the variance in g_{s} under diverse environments, then its two parameters (or three, when generalized to handle water stress) are likely candidates. If we should find that even these parameters change with lifetime or seasonal development, then we should look for a very few fundamental parameters that describe the development of the BallBerry parameters.
Critical elements of a model.
This brings up a large topic in modelling that is underappreciated, but which I
hope I can present briefly and usefully here.
When we set up a description of phenomena, there are several key
elements, some of which seem similar but with distinctions that are important
and useful. One element is the set of response variables  also called state
variables. For a model of phenomena that
occur over time, such as plant growth, these are the variables that evolve in
time. Other items also vary in time,
such as weather, but we do not predict their development. We see them as prescribed by things outside
our system, and as unaffected by the responses made by our system. We call these driving or forcing variables  e.g., air temperature, windspeed,
etc. as functions of time. Another
element of a model is initial conditions
 in what state did the system begin?
This can have major effects on the outcome, so that we need to specify
them with some care. Yet another element
is the set of parameters. These are, as noted earlier, the unchanging
elements in the processes. The example
above suggested that, at least in some cases, the two parameters of the
BallBerry equation are parameters.
Another element is a set of boundary
conditions  relations among the state variables that are maintained at the
outer limits of our system. A clear
example is air temperature far above or upwind of the canopy. Air temperature within the canopy is a state
variable  the energy balance of all the leaves heats or cools the air, in a
way that varies at all times. However,
far from the canopy (and at a fixed reference point), air T must reach the
single value measured at the weather station.
Finally, or, really, right from the start, there are the process equations  the equations that
describe each and every process that affects the state variables. In the models we will use here, these include
the enzymekinetic equation of photosynthesis (the Farquhar  von Caemmerer 
Complex models vs. simple models. This should not be "versus," really. The models can be used in fashions that are complementary to each other. I describe in this document two major models of canopy photosynthesis and transpiration, and their responses to stress. The complex one uses a great deal of detail about the structure of the canopy, while the simple one models only the sunlit and shaded leaf areas as uniform entities. The relation of these sunlit and shaded areas to detailed canopy structure is set, for one particular canopy structure. This done, the simple model runs much faster and is easier to comprehend. However, it cannot be applied with high accuracy to new canopies of different structure, unless one runs the complex model at least once to parametrize the simple model again. This need is a manifestation of the need for complex models if one is to use them in arbitrary conditions, or, we may say, to have the models be transferrable between sites and conditions. The extreme case of nontransferrability is the use of a purely statistical model, a fit to data that applies to one site and one limited set of conditions. We cannot do this here.
The complex model is "datahungry," requiring much more information to use it. This may be a realistic expectation  canopies (or systems in general) differ in many details. Some of the details are important for the results that a user is focusing upon, others are not. This leads to another use of complex models  determing which descriptors of the system are important to the results (simulations, predictions) being examined. One can run the complex model with variations in each descriptor, say, foliage density, or rootlength density, or average air temperature, and see how much difference each factor makes. For the factors that don't matter much, we can set them as constants in a simple model or otherwise make them unnecessary to specify.
There remains a hazard in complex models, that of compensating errors. A complex model may describe very many processes, each with descriptions (such as root length density) that may be hard to get from experimental data with a level of effort that is affordable. One may make guesses for the poorlyknown descriptors, and possibly "tweak" them all to get the right results for a small set of final variables. The results may have come out well only because errors in one description cancelled those in another (or several others). The only way to check for full consistency is to get a wider array of results  say, not just total growth or total nut yield and total water use, but many details of the time courses of transpiration, etc., or more deeply yet, the responses of various leaves. If these data are not obtainable with the effort that one can mount, then then one must live with reservations about the full validity of the model. We all cope with such reservations of all our models, whether they are encoded in a formal mathematical model or are nebulous and implicit assumptions we make about how the system works. In any case, the solution to poor predictions is almost never to add more processes and complexity. Judgment is needed.
Example of final use  a first estimate of optimal deficit irrigation.
Getting back to the ultimate goal, the identification of an optimal irrigation plan: I can give a preview of how the results of a model can be used to do this, mathematically. It is not a definitive model, which requires all the upcoming effort to be completed, but it shows the philosophy behind the ultimate optimization strategy.
I have made a somewhat simplified model of all the processes in a pecan orchard. This model (the "bigleaf" model described later) was constructed so that I could force the transpiration rate to be onehalf the rate of a fullywatered orchard. In particular, this was achieved by finding, by simple search, the magnitude of the BallBerry slope parameter, m_{BB}, that would do this. This choice of m_{BB} then drove a specific change in photosynthetic rate, and, consequently, in nut yield. I also resolved two stages in the growing season  leafing out, and nutfill. In the leafout stage, reduced water availability also reduces the development of leaf area, thus, of orchard photosynthetic capacity. I ran the model not only for one size of the reduction in transpiration (50%), but a whole range of the transpiration fraction, f. Furthermore, I considered different levels of reduction in the two parts of the season. After all, 50% reduction in both parts might not be best. Perhaps it is best to cut more severely during leafout, getting a smaller canopy, but supporting that canopy more fully as it is used for nutfill.
Here are the premises of the simple
optimization model:
* Leafout and nutfill are two distinct periods, with fixed durations  call them n_{1} and n_{2}, measured in months.
* Only photosynthesis (PS) during nutfill contributes to nut yield, and nutfill is a fixed proportion of total PS during this interval. That is, there's no significant reallocation of photosynthate with different stress levels.
* In leafout, stress as a reduction in wholecanopy transpirational rate (to a fraction f_{1} of normal) decreases the development of total leaf area, as the leading effect.
* In nutfill, stress as a reduction in transpirational rate per unit leaf area (to a fraction f_{2} of normal) decreases photosynthesis per unit leaf area, as the leading effect.
* Water use rates per unit leaf area (E_{1} and E_{2}) differ between the two intervals. Weather during leafout is hotter and drier, so that E_{1} is higher than E_{2}.
From these premises, we can compute
estimates of:
* Total water use, as a function of the chosen irrigation fractions f_{1} and f_{2}.
* Total photosynthesis, hence, total nut yield, as a function of f_{1} and f_{2}.
Finally,
we can find the optimal strategy of deficit irrigation:
* We can normalize these by the total water use with full irrigation. Call this F_{W}, for fraction of normal water use.
* We can also normalize nut yield the same way, dividing by nut yield under full irrigation. Call this normalized yield F_{Y}.
* Let's set F_{W} to a desired (or, really, required) value, say, ½ .
* Then, for any choice of the fractional irrigation rate in interval 1, f_{1}, there is a matching value of f_{2} that makes F_{W} = 0.5.
* For these values of f_{1} and f_{2}, we can compute (from the earlier premises) the value of F_{Y}.
* Finally, we can find the best value of f_{1} (and, hence, also f_{2}), as the one that makes relative yield, F_{Y}, the highest. With some of the simpler formulas I give for all the processes above, one can get an algebraic equation that is directly solvable. In general, we can find the optimum numerically, by simply varying f_{1} in, say, a spreadsheet (which I provide, shortly), or a Fortran program, and seeing for which value we get the highest F_{Y}.
We will now make an estimate,
putting some quantitative expressions into the premises:
* Leafout lasts two months (n_{1} = 2), nutfill lasts 4 months (n_{2} = 4). However, we want to use rates of water use (E) and photosynthesis from the peaks of each period. In leafout, water use is less than that for 2 months at full cover, because leaf area is not at final value all the time. Let's make an effective value of n_{1} that is lower, say, n_{1,eff} = 1.5. Similarly, the later part of nutfill is cooler, and leaf nitrogen is declining. Both E and PS are cut. Let's use an effective value of n_{2}, similarly. I propose n_{2,eff} = 3.
* Reduction of leaf area (cover) by deficit irrigation during leafout: While I have no direct data, I made a guess that a 50% cut results in a 30% drop in final leaf area (less than proportional  same idea as the lessthanproportional drop in photosynthesis). Let's make a linear estimate, that the fractional cover (fractional leaf area development), f_{cover}, responds linearly to the irrigation fraction, f_{1}:
_{}
We can easily set a and b from f_{cover} reaching 1 when f_{1} = 1 (so, a+b=1) and from f_{cover} reaching 0.7 when f_{1} = 0.5 (that is, 0.7 = a*0.5+b = a*0.5+[1a]). The last equation gives us a = 0.6 and b = 0.4.
* Increase in wateruse efficiency under stress, in the nutfill interval: First, let's define WUE as the ratio of photosynthetic rate to transpiration rate. The results of my more detailed simulations show a pronounced rise, hitting the limit of 1.8 times higher than "normal" under zero irrigation. Taking the response, again, as nearly linear, and denoting the value of WUE under full irrigation as WUE_{0}, we get
_{}
This form gives us the proper ratio at f_{2} = 1 and f_{2} = 0.
* To get the photosynthetic rate (per leaf area) at reduced irrigation, we multiply transpiration rate (also per leaf area) by wateruse efficiency. We have defined the stress factor f_{2} as the relative rate of transpiration per leaf area. This gives us the photosynthetic rate per unit leaf area, A, as
_{}_{}
* To get the wholecanopy rate, we need to multiply by the canopy leaf area, which has been reduced by deficit irrigation during leafout. Thus, this rate is equal to the normal rate, multiplied by the last two factors in the previous equation and by the factor (0.6*f_{1}+0.4) from the earlier note. The complete factor of proportion is then
_{}
* This captures the effect of deficitirrigation fractions, f_{1} and f_{2}, on total yield. We now need to figure out how these irrigation fractions affect total water use, relative to full irrigation.
Total water use in leafout is the normal rate of water use (call it WR_{1}), multiplied by the relative reduction in canopy transpiration, multiplied by the effective duration of leafout. The relative reduction in wholecanopy transpiration is the relative reduction in rate per leaf area (simply f_{1}) multiplied by the relative reduction in canopy area:
_{}
We can do the same for the nutfill interval:
_{}
The total use in both periods is the sum of these. Using the estimate that WR_{1} is 1.5 times larger than WR_{2}, we can write the total water use, W, as
_{}
The rate when full irrigation is applied is found by plugging in f_{1} = 1 = f_{2}. Call this W_{0}:
_{}
The relative reduction in water is the total use with deficit irrigation, divided by the total use with full irrigation:
_{}
The factors of WR_{1} cancelled out between numerator and denominator.
* As noted earlier, we are given the required cut in F_{W}  say, 0.5. Now we can choose f_{1} and see what f_{2} has to be. For example, if f_{1}=0.5, then f_{2} has to satisfy
_{}
We can do the reverse, too  we can specify f_{2} and ask what f_{1} has to be to reduce total water use by half.
* The final task is finding, for any choice of f_{1} and thus of f_{2} (or vice versa), what is the effect on yield, as the value of F_{Y}…and then testing various values of the f's to get the highest F_{Y}…for a specified cut in total water use. Here is a table from the spreadsheet I made:
f2 
f1 
fcover 
WUE2/WUE2un 
F1 
F2 
FY 
1 
0.434701 
0.66082 
1 
0.287259 
0.66082 
0.66082 
0.9 
0.48732 
0.692392 
1.08 
0.337417 
0.623153 
0.673005 
0.8 
0.542961 
0.725776 
1.16 
0.394068 
0.580621 
0.673521 
0.7 
0.601695 
0.761017 
1.24 
0.4579 
0.532712 
0.660563 
0.6 
0.663575 
0.798145 
1.32 
0.529629 
0.478887 
0.632131 
0.5 
0.728633 
0.83718 
1.4 
0.609997 
0.41859 
0.586026 







0.85 
0.514758 
0.708855 
1.12 
0.364889 
0.602527 
0.67483 
In the last line, I made a rough optimization. I saw that the best F_{Y} occurred between f_{2} = 0.9 and f_{2} = 0.8. At the intermediate value of f_{2} = 0.85, we require f_{1} to be 0.515. Now, note that this means applying 51.5% of the water that would be transpired by the smaller canopy, with only 70.8% of the normal area. Thus, during leafout, we apply 0.515*0.708 = 36.5% of normal water amounts per ground area (that's the figure in the column labelled F1 in the table). We still get 68% of normal yield. Note that f_{2} = 0.85 means, similarly, applying 85% of what the reduced canopy area would use at full tilt. The leaf area has been reduced by the factor f_{cover} = 0.708. Thus, the water application rate is 0.85*0.708 = 60.2%
If we chose to apply 100% of normal water during leafout (f_{1} = 1), then we'd have to cut F_{2} to 12.6%! (Calculations not shown.) This would give us only 21% of normal yield.
I hope that this introduction to the specific problem, to modelling, and to generating an optimal strategy shows you the importance and the promise of the whole research effort. I hope also that it sets the stage for seeing
* what specific stress responses we must characterize (canopy development, photosynthetic capacity, stomatal control, …),
* how carefully we must design the experiments,
* and how carefully we must execute each measurement.
After all, this is a simplified model. We have to verify that the parameters are right, in detail, and we have to use better representations of changes in cover, E, and A than the simple linear responses used here. Still, we see that the "obvious" strategy (water fully during leafout and cut back later) is highly unlikely to be the best.
1) Develop models that are specifically
designed for the purposes noted: synthesizing our understanding of the
stress responses; helping design experiments; and contributing to creation of a
decisionsupport system for deficit irrigation
2) Verify models with new and past field
measurements / design the best experimental protocols
3) Use the models to parametrize an initial optimization model; specifically, generate parametric descriptions of the waterproduction function of pecan orchards and of how their wateruse efficiency varies with stress level. Several auxiliary tasks include generating descriptions of Nstress responses.
Details of Objective
(1): Develop the models
Model canopy photosynthetic rates (A), transpiration rate (T), soil evaporation rate (E), and , canopy N content at fixed seasons, with constant (average) weather.
The need for
both complex and simple models
In view of all the different stress responses that we are aware of, we might start with a very comprehensive model.
In view of the need to get to a usable optimization model, and a usable decision support system for growers, we need to collapse detail to a much simpler model. This means taking the results of the complex model and fitting the simpler model to them. We'll get to the details later.
My former graduate students and I have developed a complex model. It is coded in Fortran 77 (E7E.aragorn.f), and there are over 1000 lines of comments that explain what it is doing  what the processes are, and how the mathematical solutions are obtained. Here I give a narrative version, less quantitative. You can see the quantitative information in the code comments.
The importance of canopy structure. The model resolves the actual structure of the orchard, in which leaves get different light levels, because: 1) they present themselves at different angles relative to the direction of the sun, and 2) other leaves intercept part (or all) of the direct sun and part of the diffuse skylight; this includes leaves on neighboring trees, a complexity first addressed in 1976 by Norman and Welles. We generalized it, so that all trees can differ in location of their centers and in height and shape (modelled as ellipsoids). I also used a very accurate representation of how light reaches leaves. The direct solar beam arrives statistically at any spatial location, with a probability calculated by Beers' law and using the real distribution of leaf angles and the total possible obstruction by leaves on all trees between the sun and the location. That is, there is a probability P_{dir} that the direct beam arrives at full intensity and a probability 1P_{dir} that is completely blocked at this location. The diffuse beam arrives determinstically, at a fraction P_{diff} that is also computed from Beers' law, but applied to beams from 25 different sky directions. The essential determinism is a fact I found in modelling in 1984.
At each location, there are leaves with different angles relative to the direct solar beam. We resolve 10 ranges of angles, thus, 10 ranges of direct solar radiation.
We sample the environment of leaves at 125 (or, an adjustable number of) locations within the canopy. Each location is representative of the same volume of canopy (same number of leaves, and same leaf area) as every other location, by a judicious choice of geometries.
Photosynthesis by the whole canopy. The total photosynthetic rate of the canopy is the sum of the rates for all the leaf areas. It would be computationally very inefficient to compute separately the rate for each location and each leaf angle and each class (directly lit or not). Instead, we add up, over all locations, the fraction of leaf area in 10 ranges of total light level (called irradiance). Then, we compute the photosynthetic rate (and transpiration rate) for the 10 different irradiances (the midpoint of each irradiance "bin"). Now, leaves at different locations also see different temperatures, windspeeds, and humidities. We ignore this as a simplification. There are possible ways around this, but they are very demanding, not only in computation, but in describing how heat and humidity and wind move around in the canopy.
The complete environment of the leaf determines photosynthesis. The leaf photosynthetic rate, A_{leaf}, depends not only on the irradiance (in photosynthetically active radiation between 400 and 700 nm in wavelength), but also on temperature, humidity, CO_{2} concentration, and windspeed. There are 3 (really, 4) basic equations that capture the greatest part of the biophysical and biochemical responses and allow a computation of A_{leaf}….and, at the same time, the leaf transpiration rate, E_{leaf}, and the stomatal conductance, g_{s}:
1) The Farquhar  von Caemmerer  Berry model of A_{leaf} in terms of basic photosynthetic capacity (V_{c,max}, related to content of Rubisco enzyme, in essence, and closely related to leaf N content), CO_{2} partial pressure at the Rubisco sites (C_{i}), and leaf temperature.
_{}
Here, Γ and K_{CO} are temperaturedependent functions for the Rubisco enzyme.
To use this equation, we need the latter two quantities. We get these from the next two equations:
2) The equation of energy balance. We sum up all the ways that a leaf can gain and lose energy, and assume that the leaf is always close to steady state: the sum is zero. On average, this is an excellent approximation.
We account for energy gain from radiation  the PAR portion of the spectrum (close to half of solar radiation), the nearinfrared portion (NIR; the other half of sunlight), and thermal radiation. We have already computed how much PAR reaches various amounts of leaf area. We approximate that the same amount of NIR energy reaches these leaves. This is a weak approximation, because NIR is absorbed much less strongly; it bounces around in the canopy and reaches leaves deeper in the canopy. This bouncing also means that a significant amount of NIR reaches leaves after first scattering off other leaves.
The thermal infrared radiation (TIR) arrives from two main sources the sky, radiating from water vapor molecules at a range of altitudes, and the other leaves. We ignore the soil, for now. TIR is radiated by a "perfect black body" at a rate proportional to the fourth power of the absolute temperature of an object. A real body has this rate, multiplied by an emissivity (nearly 1 for leaves). The sky has an effective temperature that is a convolution over temperatures of the air at all altitudes with the amount of water vapor at each altitude. We can make a good approximation that it has an effective temperature that is a fixed number of degrees below air temperature at a weather station. In the code, we assume that all leaves have the same temperature, the average in the canopy. The sky TIR radiation penetrates to a given location in the same proportion as diffuse sunlight.
We also account for evaporative cooling of the leaf by its own transpiration rate. This depends on the stomatal conductance (and a larger boundarylayer conductance, in series), the leaf temperature, and the partial pressure of water vapor in the surrounding air. Now, the leaf temperature is part of what we're solving for, so we must solve the energybalance equation iteratively. The boundarylayer conductance depends, by a simple formula derived by others, on the leaf linear dimension and on the windspeed at its location. We assume that it is the average windspeed (really, about ½ of this would be better). The stomatal conductance must be estimated from yet another equation:
3) The BallBerry equation for stomatal conductance. This was also mentioned above. It takes the simple form
_{}
Here, A is the leaf PS rate, A_{leaf}, and h_{s} and C_{s} are the relative humidity and the CO_{2} mixing ratio at the leaf surface, beneath the leaf boundary layer. The occurrence of A in this equation means that this equation must be solved iteratively with the photosynthesis equation. First, this represents real physiological feedbacks and feedforwards, as noted earlier. Second, there are efficient ways to do this solution, which I can describe shortly. The values of h_{s} and C_{s} depend on A and E of the leaf and on the stomatal and boundarylayer conductances, using readilyderived formulas of physical transport phenomena.
4) The equations of transport for CO_{2} and water vapor: there must be a drop in CO_{2} partial pressure from outside to inside the leaf, to support the inflow of CO_{2} during photosynthesis. There must be a rise in water vapor, to support the transpiration rate. This gives us the relation for internal CO_{2} (C_{i}, the quantity need by the PS enzymekinetic equation) to external CO_{2} partial pressure, C_{a}:
_{}
Here, P_{air} is the total air pressure (this comes in because we use socalled molar units of conductance, the modern convention, for various good reasons), and g_{tot,CO2} is the total conductance for CO_{2} through the stomata and the boundary layer.
Solving all these coupled equations. The solution method I developed is efficient, using a binary search over magnitudes of g_{s} until all the equations are solved simultaneously. Consider a guess at the value for g_{s}. The energybalance equation has all the other quantities specified, so this lets one calculate the leaf temperature (with quick iterations). We combine the enzymekinetic equation (1), with its parameters corrected for the leaf temperature, and the transport equation (4) to get a single equation for C_{i}. When we use the form of the enzymekinetic equation generalized to handle lightlimited photosynthesis, this becomes a quartic equation. We solve it rapidly by iteration. Now we have both C_{i} and A. Finally, we rewrite the BallBerry equation to highlight the error in the solution, as
_{}
When we have the right guess for g_{s}, F becomes zero. We home in on the proper value of g_{s} by a binary search. We guess the min and max values that g_{s} could lie between. We compute F at each end, and then for g_{s} in the middle. The solution has to lie between the values of g_{s} where F changes sign. We take these two values as the new min and max, thereby halving the interval. We keep doing this until the interval is less than some preset accuracy, say, 0.00001 mol m2 s^{1}. Only about 15 iterations (halvings of the original interval) are needed for this.
Accounting for respiration. A final note about photosynthesis is that the gross rate is debited for instantaneous respiration in the leaf. This has been found repeatedly, including by us, to be 8 to 10% of gross PS at the current twoweekaverage air temperature, T_{mean}. We input the latter and calculate the respiration rate for any leaf, applying an exponential factor in actual leaf temperature, exp(0.07*(TT_{mean})).
The rate of photosynthesis is not to be compared with net CO_{2} exchange of an orchard, because respiratory losses of CO_{2} (partial undoing of photosynthesis) occur at night everywhere, and at all times in the trunk and in the soil…at a rate that makes net CO_{2} uptake as small as 20%, or even 0% or less, of this "canopy gross" photosynthetic rate. The soil respiration is typically largest. It comes from living root tissue, when sugars are metabolized for energy to drive synthesis of new tissue and to maintain all tissue. It also comes from microbes in the soil, using up direct exudation of sugars and acids by the roots (done by the tree for a variety of ecological reasons) and also breaking down dead roots, which arise on a short turnover time from live roots.
The canopy affects its own environment. We're not done. The transpiration by all those leaves adds humidity to the canopy, changing the environment of the leaves. So, too, does photosynthesis drop the CO_{2} level, and convective energy transfer alters the air temperature in the canopy. We have to iterate the solution for, particularly, the air temperature, T_{air}, and watervapor partial pressure, e_{air}. At each iteration, we get a new e_{air} and a new T_{air}…and then new canopy totals of A and E…which gives us new e_{air} and T_{air}. The iterations are prone to oscillate and diverge, so that I apply a limiting factor. The changes in e_{air} and T_{air}, from their values in "free" air above the canopy, depend on the boundarylayer (or aerodynamic) conductance of the canopy as a whole. This depends inversely on windspeed, with a constant of proportionality that depends on canopy leafarea index, LAI. We know windspeed from the weather data, and LAI by any of various measurements or even some robust rules of thumb.
All these processes change in rate over the day, as the solar angles change and so does air temperature, humidity, and windspeed. We recompute all the quantities above at any time step desired (currently, 10 minutes, to match some sapflow measurements we made in the field. We are changing this to 1 hour, to match weather data and reduce "overkill" in detail).
Some processes are still left out. With this complexity, we still have left out a number of processes:
* We don't look at the energy balance at the soil and its related soil evaporation rate (but we are adding this, as we are also adding it to the next simpler model, described shortly).
* We assume a constant canopy photosynthetic capacity, while a real canopy declines in leaf N and this PS capacity as N is withdrawn to support protein synthesis in developing nuts (but we are adding this also).
* We take the canopy as dry, not accounting for interception of (very modest) rainfall and its subsequent evaporation; we have no current intention of adding this process.
* We take the stomatal control parameters, m and b, as constant. Really, under water stress, m certainly declines, and I modelled this in other work, using field data. We are putting in this process, in two separate versions of the complex model. In one version, we use the real soil water balance, computed using our estimates of wholetree transpiration and the inputs from rainfall and irrigation. This is used to compute soil water potential. It is also used to compute the hydraulic resistance from bulk soil to roots; this, multiplied by wholetree transpiration, gives us the drop in water potential at the roots. The root water potential may be used to estimate the drop in BallBerry slope, m. In a second version of the model, we simply enforce a halving (or other change) in wholetree transpiration. We search for the value of m that achieves this. There are consequent changes in leaf temperature, g_{s}, and A, all of which we use to see how wholetree photosynthesis is changed when water use is halved.
The second version of the model leads fairly directly to the optimized irrigation strategy. The first, more complicated version, resolves what happens in real time with real weather and irrigation. Thus, it is useful in a decisionsupport system (although one would not want to run it in real time for any grower's conditions; we will use it to derive simpler models). Also, the first version resolved details of photosynthetic rates, for example. As such, it generates testable hypotheses about how exactly the coupled processes all change under stress. Field measurements will be used for this testing, which is critical to affirm our understanding…or correct it.
A big leaf. Much of the complexity in the model above arises from describing the details of how sunlight penetrates to leaves in all different locations and in all different orientations. Many researchers have used a surprisingly good simplification, called the bigleaf model. It treats the leaves as falling into two classes, sunlit and shaded, all of them in each class having the same irradiance and other environmental conditions. For a given leafarea index, LAI, one can estimate from very good models of light penetration (in uniform canopies, usually) how much leaf area is directly sunlit, at a given solar elevation. For example, with overhead sun and a very deep canopy of randomlyoriented leaves, the sunlit LAI is 2. There are ways to approximate the variation with finite total LAI and for different solar elevations (I need to do the latter better).
Solving the same equations. For each class, sunlit and shaded, one computes A, g_{s}, and E per unit leaf area using the same set of 4 simultaneous equations as used in the complex model  that is, the enzymekinetic equation for photosynthesis, energy balance, stomatal control, and transport. For transport, one has to estimate an average total conductance from leaves to free air, from the canopy aerodynamic resistance and the boundary layer resistance of individual leaves.
Programming, and accounting for nutfill. I have programmed this model in Fortran 77, as the file ETstressed_full2.f. It already includes the effect of nutfill drawing N out of the leaves and reducing photosynthetic capacity. It uses the N demand in nut kernels themselves, inflated by 50% to account for growth of the sheath and the petiole. It assumes, as Ted requested, that no new N is taken up from soil, so that all this N is drawn from the leaves. The PS capacity, V_{c,max}^{25}, is reduced in direction proportion to the N fraction in the leaves. (The capacity is referred to a common reference temperature of 25ºC; at any chosen single temperature, it is proportional to the N content.)
The simple model also converts the gross magnitude of wholetree photosynthesis to the rate of growth of nuts. It assumes that a constant fraction of photosynthate (newlymade sugars) ends up as nut biomass, after losses for biosynthesis, maintenance, and turnover to tissue. The fraction has been adjusted to replicate realistic nut yields. Its value is also consistent with known patterns of biosynthesis, maintenance, and turnover, and with known fractional allocation of growth to nutfill.
Using only one time of day as representative. The model computes PS and transpiration for a single time of day, typically midday, and then scales the rates to the whole day. One multiplies the rate by the length of the photoperiod, naturally. That is not sufficient, because both rates are lower away from midday, as sunlight interception is lower then, and so is temperature. A simple estimate is that the rate varies as a sinusoid from dawn to dusk. A better estimate takes into account the real variations of light interception and temperature, seen in the complex model.
This simplification of the diurnal course of canopy fluxes and nutfill can readily be relaxed. The model can also be reprogrammed with little effort to perform simulations over a full diurnal cycle.
The model allows the user to enforce any desired cutback in water use  the factor f (f_{1} or f_{2}) noted in the description of optimal deficit irrigation, above. It does so, again as noted earlier, by searching for the value of the BallBerry slope that forces transpiration to be half of the unstressed (fullywatered) rate. The consequences for photosynthesis and nutfill follow from the effect of the new stomatal conductance on all the coupled processes. I should note in this regard that the reference rate is not corrected for the change in canopy leafarea development from prior stress. This will have to be incorporated in a better version of the model.
Using the models, "forward" and "inverse"modes
Let's start with model role (2) described earlier  designing experiments. A major part of this effort is to: 1) make our best initial guess of how the plants respond to stress, and the rest of the environment; 2) run the model (either one, complex or simple) with this guess; and 3) see what effects of stress are most important and most pronounced.
Here is a logical set of steps:
* Obtain the inputs needed by the model. Understand what they signify, and know the units of each input. Find the best data for the purpose.
* Run the model (learn the mechanics of writing the input file, executing the code)
* Examine the outputs. Understand the items that have been computed, including their units, which give their spatial and temporal scales, among other things. Address several key issues:
 What patterns in the variables, temporally, spatially within the canopy, and as sensitivity to inputs, particularly to physiological traits and meteorological conditions? This type of analysis is needed in using models to design both the field experiments and the decisionsupport systen.
 How does one distill the results to simpler quantities that we want so that experiments can test the accuracy of our understanding (our model, as the embodiment of all we know about the system). For an example, play with the simulation "data" to determine the most productive way to characterize the fractional reduction in leaf PS at various stress levels.
 How do we find out which physiological and treestructural parameters are most important to measure accurately? That is, which ones show the greatest effect on nut yield and water use if they are varied within the bounds that are typical?
We will now consider each of the steps in some more detail.
Inputs to the complex model: sources of
data, quality control, formatting, and units
A note on potential changes in required inputs. The model will be revised during the current work, to describe water stress and to make it more userfriendly for the new users.
Meteorological and other environmental information. Here are the entries, in the order in which they are input to the program (though there are other data items that precede them or that intervene between some of these items):
* weafile the name of the file that has the hourly weather data that covers the period of the simulation.
Format: text (characters), up to 48 characters long. No spaces can occur in the name (pet peeve: Apple started the practice of allowing spaces in file names, with its very abundant pitfalls in computing).
The contents of the file are editeddown entries from the standard data on the site, weather.nmsu.edu.This file should be in a fixed format  that is, the exact spacing of entries is observed. I can modify the program to allow freer formatting, just as long as the entities appear in the correct order; the fixed format was easily created by standard editing on a UNIX machine. Each line in the file covers one hour. Here are the first three lines of a file that I used to represent weather at the Hort Farm in July, 2003 (they are not indented in the original):
07 1 03 01AM 74.6 55 00.00 1.0 000 00.00
07 1 03 02AM 72.3 56 00.00 1.0 000 00.00
07 1 03 03AM 71.2 56 00.00 1.0 000 00.00
The entries, all but two separated by spaces, are in the same order as on the weather site. They have been pruned down, retaining only the following:
month, as a 2digit integer  with a leading zero if it is less than 10
day, as a 2digit integer
hour, as a 2digit integer  again with a leading zero if it is less than 10
ampm  indicator of AM or PM; note that 12 AM actually is noon, contrary to standard use
air temperature, in degrees Fahrenheit, as a real number with 3 digits before the decimal point and one digit after  converted internally to SI units of degrees Celsius
relative humidity, in percent, as a twodigit integer
windspeed, in miles per hour, as a real number with 2 lead digits and one digit after the decimal point  converted internally to SI units of meters per second
(Not used: wind direction, in degrees clockwise from N, as a 3digit integer)
insolation (solar energy falling on a horizontal surface), in MJ per m^{2} per h, as a real number with two leading digits and 2 digits after the decimal point
* date0wea  the date on which the weatherdata file begins, which is not necessarily the date on which the simulation begins (the user may have chosen a convenient wholemonth block for the weather file, for example)  as text between quotes, such as "070103". Hyphens are required between the parts of the date, which must be in monthdayyear order. Each element must be a 2digit integer  e.g., 07 for July, not simply 7.
* Ca  the partial pressure of CO_{2} in the ambient air, in
Pascals. It affects photosynthetic rate,
hence, also, stomatal conductance and transpiration. This is not measured at
most weather stations, so you can compute it from total air pressure and an
assumed fixed mixing ratio of CO_{2} in air (today, about 385 ppm; it
varies notably only during inversions).
Example: the air pressure in
* Pair  total air pressure, in Pa; see Ca, above
Physiology and other leaf traits. Also in order of appearance, but with other items preceding the first one and items intervening among these:
* Vcmax25  photosynthetic capacity of leaves per
unit area, taken as uniform throughout the canopy, in units of micromol
(μmol) CO_{2} fixed per m^{2} per s, the same units
reported by all common measurement systems, such as the LICOR 6400. This is the rate that photosynthesis would
attain at high (nonlimiting) light intensity and high (nonlimiting) CO_{2}. It is used in computing the PS rate at any
light and CO_{2} levels, as seen in earlier equations. It is in direct proportion to the amount of
Rubisco enzyme per unit leaf area, Rubisco being the ratelimiting enzyme at
high light levels. We have some earlier
measurements of this capacity in pecan leaves, obtained by my own group in a
Salopek orchard and more extensively by David Johnson working at the
Leyendecker farm and on the new trees on
_{}
Leaf temperature, T_{L}, in ºC, is used to compute the two enzymekinetic parameters, Γ and K_{CO}, via standard formulas. This gets a little complicated, so that I have provided a spreadsheet, computing_PS_capacity.xls, that my graduate student also uses. The formulas are:
_{}
_{}
The PS capacity changes with temperature. To make it convenient to use at any leaf temperature, it is entered as its value at a standard reference temperature of 25 ºC and then its value at any other temperature is computed from this:
_{}
The value of Vcmax25
will vary from leaf to leaf, so we need to take a judicious average. Its average over leaves that get mostly full
sun is the important quantity. I propose
that we measure Vcmax25 for different
leaves that experience a variety of different lightlevels in the long term 
say, from different positions in the canopy.
Ülo Niinemets in
* thetaPS  curvature in the transition between lightlimited and lightsaturated PS rates. At low light levels, the PS rate, A, is proportional to the light intensity, I_{L} (irradiance on the leaves, in the PAR region of the spectrum):
_{}
Here, the subscript "LL" denotes "lightlimited." At high light levels, PS is lightsaturated, at the rate expressed by the Farquhar  von Caemmerer  Berry rate, above. The transition is not sharp, that A is at the lightlimited rate until it reaches the saturated rate and then suddenly is constant. A very good representation of the transition is
_{}
Here, A_{max} is the lightsaturated rate.
The value of θ can be derived by making measurements of the photosynthetic rate at a variety of light levels on single leaves and then doing curvefitting This can be tedious. Among many species, θ takes values near 0.8, so we use this.
* RdovrA  the ratio of leaf respiration to photosynthesis in the short term, as a dimensionless number such as 0.08. Respiration occurs all the time, in leaves, stems or trunks, roots, and soil. This factor accounts for the immediate daytime respiration. The PS rate that we measure with a gasexchange system is the net rate, after respiration. We formulate the net rate of PS, A = A_{net} as the gross rate minus respiration, A_{gross}  R_{d}, and then as A_{gross}*(1R_{d}/A). As noted earlier, it has been found repeatedly that RdovrA is in the range 0.08 to 0.1. I use 0.08 as a good estimate. (If you wish, you can measure this in the field, too.) This is the ratio when leaves are at the longterm (2week) average of the temperature over the photoperiod. That is, leaf respiration acclimates, not rising with temperature in the long term (or, not any more than does photosynthesis). In the short term, during a single day, respiration does vary with temperature in an exponential fashion,
_{}
* aPAR  the absorptivity of leaves for photosynthetically active radiation, as a dimensionless number, say, 0.85. This can be measured in the field with a variety of instruments. It does vary with leaf N content. The value of 0.85 is a good average value.
* mBB  the slope of the BallBerry relation for stomatal conductance, as a dimensionless number. This is constant over the medium term, but it may vary seasonally; we need some gasexchange measurements over the seasons to determine this. A mean value is about 10. The value has a significant effect on PS rate and wateruse efficiency.
The value of mBB can be extracted from a series of gasexchange measurements on leaves in different environmental conditions (light levels, temperatures, humidities). This series may be run on a single leaf, kept in the chamber while the environment is changed, using controls of the LICOR 6400. It can also be done from measurements on different leaves found in different conditions at any one time of day. One has to assume that the leaves share the same value of mBB (which is an assumption made in the full model, anyway.
One computes the BallBerry index, A h_{s}/C_{s}, from the measurements at each condition, and then plots the values against the values of stomatal conductance, g_{s}, at each condition. The value of A comes directly from the LICOR. The relative humidity at the leaf surface beneath the boundary layer, h_{s}, is offset from the value in free air. It involves the boundarylayer (g_{bL}) and stomatal conductances and the watervapor partial pressures in free air (e_{air}) and in the leaf (e_{L}, computed from leaf temperature). The equation is
_{}
Another equation relates C_{s}, the mixing ratio for CO_{2} in air at the leaf surface, to the mixing ratio in free air, C_{a} (remember, not in Pascals), the photosynthetic rate, and the boundarylayer conductance:
_{}
These calculations are not programmed into the LICOR. One must import the data from the LICOR into a spreadsheet. I am supplying as a separate file a spreadsheet, BB.xls, that makes all these calculations from data that are readily pasted in from a set of LICOR runs. It also calculates Vcmax25. It uses data that David Johnson acquired in his M.S. research in 2003. The many columns of LICOR data will become familiar to the user, with practice. I don't have the time or space to explain all of them here.
* bBB  the intercept
in the BallBerry relation, in mol m2 s1.
This is obtained in the same datafitting exercise as is mBB, above. In the absence of new data, a value of 0.02
is a good estimate.
* dleaf  the characteristic dimension of the leaf, usually the short side, in meters. It is used in calculating the boundarylayer conductance of the leaf as 0.264*sqrt(u/dleaf), where u is the windspeed in m/s. A simple visual estimate will be sufficient. I used 0.05 m = 5 cm.
Canopy structure. In order of appearance, these are:
* fd  the average
foliage density in the canopy, in m^{2} per m^{3}, that is in m^{1}. This can be measured by various instruments,
including a LICOR LAI2000 canopy analyzer.
This instrument has sensors that view upward in 5 ranges of zenith
angles, at all azimuthal angles. It uses
some clever theory to relate this to the leaf area index that lies above the
position of the analyzer. At the base of
the canopy, the instrument would see neighboring trees unless one set it to
look only in a narrower "cone" upward, that is, by turning off the
sensors that look at angles closer to horizontal. We have done this analysis in fully developed
pecan canopies and have found fd =
0.6 m^{1}. Note that this is
related to average LAI over all areas of the ground. This LAI is the total leaf area in the crown
of any one representative tree, divided by the ground area dedicated to each
tree. Commonly, the trees are spaced at
30 feet apart, essentially 9 m, and if they have nearly spherical crowns that
just touch. That gives them a canopy
radius of R = 4.5 m, a canopy volume of (4/3)πR^{3} = 382 m^{3},
a total leaf area of 382*0.6 m^{2} = 229 m^{2}, and a ground
area per tree of 9*9 = 81 m^{2}.
Thus, the average LAI in an orchard is 229/81 = 2.8.
* rhoovrC1  this is a proportionality constant in computing the aerodynamic conductance of the canopy as a whole (the inverse of a resistance to transport of water vapor from the canopy to the free air). It is in units of mol m^{3}. The conductance is formulated in a welldocumented model called SiB, developed by Piers Sellers and others. The formula is g_{b,can} = rhoovrC1*u, with u the windspeed above the canopy in m/s. The value varies with the total LAI of the canopy and some details of how the leaf area is distributed by height. For a fully developed pecan canopy, a value of 1.5 is appropriate.
* ntrees  the total number of trees in the simulation, as a simple integer. The user will describe the geometry and location of all these trees, relative to a central tree whose photosynthesis and transpiration rate is being calculated. The central tree must be one of these. For each tree, the user provides the coordinates (x and y) of the center of each crown, as well as its shape. The program thereby can compute how all these neighbors interfere with light interception by the central tree. In a typical orchard, the trees are fairly uniform, so that the coordinates are simple multiples of the basic row spacing, and all the crowns are specified to have the same geometry The potential to describe nonuniform trees was important in the original use of this model at the Bosque del Apache. One does not need to specify all the trees in an orchard, which would cause the calculations of light interception to take an extremely long time. We have considered trees up to 4 rows away in both x and y (EW and NS) directions. This gives us 4 + 1 + 4 = 9 trees in a row, and 9 x 9 trees overall. This is adequate for computing light interception except at very low solar angles. These low angles occur for a very limited duration in the day, and at times when PS and transpiration are small, anyway. Consequently, they are not important.
* xtree, ytree, ztree 
the coordinates of the center of the crown of each tree, in meters. If a typical pecan tree has a crown of 4.5 m
radius that begins 1.5 m above the ground, then the center is at ztree = 4.5 + 1.5 = 6 m. The central tree is defined as being at the
origin (xtree = 0 = ytree).
If the orchard is uniform, all the other trees occur at coordinates that
are multiples of 9 m. Thus, the first
tree to the east (which is the positivex direction) has its center at (9, 0,
6). The second tree to the east is at
(18, 0, 6). The first tree to the north
is at (0, 9, 6).
* atree, btree  the
semiaxes of each tree's crown, both in meters.
These items allow the tree crown to be described more generally than as
a sphere, namely, as an ellipsoid. Its
"radius" in the horizontal
plane is atree, and its
"radius" in the vertical direction is btree. For the spherical
crowns mentioned earlier, both these quantities equal the radius of 4.5 m.
* thetac, phic  in the spirit of generalizing the geometric description of the tree crowns, these quantities describe the tilt of the crown, as a zenith angle thetac and an azimuthal angle phic. They are entered in units of degrees and are converted internally into radians. If the crown is nonuniform, and if its morevertical axis is tilting offvertical by 30º, then thetac = 30. The azimuthal angle is measured clockwise from the north. Thus, if the canopy tilts east, phic = 90.
Controls over desired accuracy and over reuse of information. Accuracy in any part of a calculation comes at a price, that of more computer time. The time for any simulation can be kept to a practical value (say, no more than 1 hour even for the most demanding simulation) with judicious choices of the error tolerances and related quantities.
* nbins  the number, as a simple integer, of discrete light levels (irradiances) that will be sampled when computing photosynthesis rates over all leaf area. Recall that the model samples many locations in the canopy and many leaf orientations. It saves time, and achieves a reasonable accuracy, by considering light levels in only a few discrete ranges, and then computing PS (and transpiration, etc.) for the average irradiance in each range. I typically use 10 such "bins". If the maximum irradiance of a leaf in full sun is 2000 μmol m^{2} s^{1}, then each bin is 200 of these units in range. The first bin goes from 0 to 200, the second from 200 to 400, and so on. Why do we need to do this, rather than calculate one grand average irradiance on all leaves and use that in calculating PS and such? It is because the photosynthetic rate is a very nonlinear function of irradiance. It begins as a straight line with a positive slope and ultimately turns into a flat horizontal line. Because of this convex shape, the PS rate calculated at the mean irradiance is higher than the mean of the PS rates at all the different irradiances encountered. Keep the value nbins = 10 for a good balance of accuracy and computing time.
* D000  the PAR flux density of diffuse skylight to the ground on a clear day, in μmol m^{2} s^{1}. The total solar radiation in the PAR band varies with solar angle but the diffuse skylight is reasonably constant. If there are few aerosols in the air, this diffuse flux density is about 200 in the stated units. You might keep this value. While it ranges from about 150 on very clear days to 400 on cloudy days, the changes do not affect canopytotal photosynthesis much.
* I00start  the maximal intensity of the
direct solar beam, in μmol m^{2} s^{1}_{. }In the clear air of
* ipre  a simple switch variable, an integer that is either 0 or 1. It indicates if the program, in a previous run, computed and stored the huge array of probabilities of light penetration to all locations in the canopy for all times of the simulation. These calculations take the bulk of the time in the model, so it is worthwhile to save them. They create a big file (one that I made had almost 112,000 lines and over 6.5 MB; this can be greatly simplified). Note that the calculations can be reused only if the simulation is rerun with exactly the same dates and times and the same canopy structure. Such reruns are actually fairly common  we may update the physiological parameters, or the stress levels, for example.
* ipensave  if the big array of lightpenetration probabilities has not been calculated yet or has not been saved from a previous run, then it will be computed anew. This variable is a simple integer. If it is set to 1, then the array will be saved.
* lpenfile  the name of the previously saved array of lightpenetration probabilities from a previous run with the same simulation times and the same canopy structure. It is simple text (a character string), up to 48 characters long, with no spaces allowed in the name, and no quotation marks around the name.
* hs0  a real number, as a fraction with magnitude between 0 and 1 (really, about 0.2 and 0.8). It has two purposes. First, if it is a negative number, then the program will use the full BallBerry equation in computing stomatal conductance, g_{s}, including the factor of surface relative humidity, hs. If it is positive, then h_{s} will replaced by hs0 in the calculation of g_{s}. This effectively shuts off the humidity response, if one wants to see how important this is. Second, the absolute value of hs0 is used as the relative humidity when calculating the leaf respiration rate at the mean temperature (see RdovrA).
* tol  the tolerance within which the value of stomatal conductance will be solved in every calculation, in mol m^{2} s^{1}. Conductance values range up to about 0.5 in these units, and I want a good estimate, so I have set this tolerance typically to 0.00001. This may seem excessive, but a binary search (see above) gets 1024fold more accurate in each 10 steps, so that setting tol so fine adds little execution time.
* chaval  a text string, up to 20 characters long, used to make changes in any of a wide variety of inputs before running the full simulation. Consider the case where you have set, say, Tmean to 28 (in units of ºC) and you want to change it to 32. You would compose chaval as a string with the first three characters as the keyword for Tmean, namely, tme, then a new value. That makes the string into tme28 (tme 28 and tme 28 and tme28. , etc. all will work, too). I have used this method of resetting in many programs. It is particularly useful if you are running the program from keyboard input, but it can be used even if all the input is saved in a file. The user can string together as many changes as desired, each on a new line. A blank line indicates the end of changes, causing the program to start doing the full simulation. For example, to change Tmean and then Vcmax25, one might enter lines (not indented) such as:
tme 28.
vcm 120.
The program will print out the keywords for the inputs that can be changed, if you type in as chaval the simple string key. It will print out the current values of all the changeable inputs if you type in just see. It will prompt you to reset every such input if you type in just all.
* stepsize  the length, in meters, of the steps taken by the program in tracing the path of the direct solar beam through the canopy when computing the lightpenetration probability. With the assumption that leaves are not clumped (a correction can be made) and that leaf orientations are random, the probability that a ray of light will penetrate a distance ds through a canopy with foliage density fd is P = exp(0.5*ds*fd). If the canopy is interrupted  say, in the space between the edge of one tree's crown and the next, we only want to count the part of the path that has foliage. With irregular canopies, the gaps can be in a variety of locations and of various lengths. To get an accurate sum of the foliage area (sum of ds*fd) along any one ray's path, we need to sample at short intervals. The error will be small if stepsize*fd is small. The total error in the probability, P, is only 1% if the error in 0.5*ds*fd is about 0.01. So, if fd is 0.6, then we should set ds to about 0.01/(0.1*0.6) = 0.17. I have used 0.1, which gives a good balance of accuracy with computation time.
* radmax  maximum length of the path to trace through the canopy, in computing the probabilities of light penetration to any point in the canopy. When we start from an arbitrary point within the crown of the central tree and trace the path of a light ray outward to find the occlusion by leaves in this crown and all its neighbors, we don't know where to stop. Perhaps there are some particularly tall trees in some directions. This case is so rare that we allow the program to stop looking for intervening foliage after moving the distance radmax along the ray. This value might be about 5 times the diameter of the average crown. For trees with 9m diameter crowns, we might set this to 45 m. (In the sample input file, I have 60; I should reduce this).
* nradius, ntheta, nphi  the number, each as a simple integer, of discrete radii, zenith angles, and azimuthal angles that will be sampled within the crown of the central tree is computing light levels, PS, and transpiration. Each location in the canopy has a different probability of the direct beam reaching it, and a different fractional penetration of diffuse skylight. We can only sample a finite number of locations. An efficient way to sample them is to use a set of different distances (radii) from the center of the crown, and, likewise, a set of different zenith and azimuthal angles. A value of 5 for all 3 quantities gives good accuracy for canopies with normal diameters and foliage densities. The program spaces out the radii nonuniformly, so that each "shell" (range of radii) has the same volume as all the other "shells."
* itoffset  a time in minutes, and a relic of the use of this program for simulating measured sapflows. The sapflows were measured every 10 minutes (the current timestep of the simulation), but the measurements did not always begin exactly on a 10minute mark (say, at 8:53 instead of 8:50 or 9:00). The offset here lets the simulations and the sapflow measurements get in sync. For our new purposes, set itoffset to zero.
* date0, datef  the beginning and ending
dates of the simulation, as text strings with monthdayyear, e.g., 070103 073103.
* minute0, minutef  the times of day, as minutes after midnight (that is, as integers), that the simulation runs each day. The program can save execution time by not simulating fluxes during the dead of night. (This would lead to some underestimates if the plants do nighttime transpiration, as occurs in some species.) If we set these to 0010 and 1430, the simulation does run during the night. If we set these to 0330 and 1170, then the simulation runs from 5:30 AM (5 hours and 30 minutes, or 330 minutes after midnight) to 7:30 PM (19 hours and 30 minutes after midnight). Note that all times are in standard time, never daylight savings  the same convention as used in the weather data.
* dtree, drow  the spacing of trees in the x and y (row and column) directions, in meters. These quantities are used to transform computed rates of PS and transpiration from the pertree values to the rates per unit ground area (say, kg water per m^{2} per hour or mm per day). For nonuniformly spaced trees, these parameters are meaningless.
Inputs to the simpler, bigleaf
model
Many of the inputs are similar. Among these are:
Meteorological and environmental conditions. In order of appearance, these are:
* Ca, Pair, weafile, Tmean  same meanings as in the complex model. Currently, however, the model uses constant weather, encoded in the subroutine getenv. It can be set up to read real data from weafile.
Physiological parameters. In order of appearance, these are:
* thetaPS, RdovrA, aPAR, dleaf  all with the same meanings as in the complex model
* Vcmax25, mBB, bBB  the meanings are the same, but these are made dynamic, varying over the simulation interval because of stress. The input values are those at the beginning of the nutfill interval.
Controls over accuracy. There is one parameter that is the same as in the complex model:
* tol  accuracy requested in solving for g_{s}.
Other inputs are novel because the specification of the simulation interval is a little different, or the canopy structure is very simplified, or the simple model actually tracks growth and nutfill:
* LAI  average effective leaf area index of the canopy as a whole. This is used to compute light interception and sunlit and shaded leaf areas, hence, photosynthesis. I set it to 4, which may be a bit high, but the results are rather insensitive to LAI above 3 or so.
* mLa  average dry mass per area in leaves, in g m^{2}. This is used to compute the total mass of N in leaves, for the accounting of N depletion during nutfill. I set it to 129 to get a proper total N content. It also is consistent with field measurements.
* betanut  this is the efficiency of converting photosynthesize sugars into final nut biomass. It accounts for the "losses" of fixed carbon along the way for respiration (generating energy for biosynthesis and maintenance), for turnover of live roots to shed roots, and for the allocation to other parts of the tree (woody aboveground parts, sheath and petiole, roots). I chose a value of 0.12. That is, only 12% of photosynthesized sugar mass gets embodied in final nuts. It matches the yield of nuts in realistic conditions.
* fNnut  the fractional N content in dry mass of nuts. This is needed to compute the amount of N that must be withdrawn from leaves in order to support nut growth. I set it at 0.02 = 2%.
* Nscale  an inflator for N demand in nutfill, to account for N used in growth of the sheath and petiole. I used 1.5. I am ready to adjust it with full data.
* fNL  the initial value of N content in leaves, as a fraction of dry mass. As this value declines with nut growth, I scale Vcmax25 down in exact proportion. This is a first guess as to the relation of N in leaves to PS capacity.
* tolE  tolerance in trying to hit exactly 50% of normal water use. This keeps the iterations to a practical level. I set it to 0.01 = 1% (I accept 49 to 51%, that is).
* JD1, JDf  Julian days of the beginning and end of the simulations. I set 182 and 271, allowing 2 months from the first of July. Real nutfill last longer but the simulation uses an effective season length, scaled from the most favorable time.
Basic structure of the code for the complex model
Many lines of comments (nonexecutable statements of information), about the purpose, structure, processes represented in the calculations, and the accounting measures such as cumulative water use
Declare variables  names, types (integer, real, character)
Declare arrays  variables with several dimensions, such as canopy locations
Set up common blocks  listing of variables that are automatically shared between routines, to avoid having to have long argument lists when subroutines are called
Set up data statement  for unchanging data, such as keywords, number of days in each month…
Begin reading of input data
Set values of physical constants
Read and process weather data  interpolating to exact intervals of the simulation, estimating derived quantitites such as PAR flux densities from total solar energy flux, and converting to metric units
Read in more data  particularly physiological parameters, and allow opportunities to change values
Calculate leaf respiration at mean temperature  first call to the routines that solve the simultaneous equations for photosynthesis and stomatal conductance
Read in data on canopy structure  and allow opportunities to change values
Set up the 125 or so locations in the central tree where leaf PS and transpiration will be sampled, and allow opportunities to change values
If not done in earlier run: compute light penetration to all sampled canopy locations
Loop over days
Loop over times of day
3 nested loops over locations
Loop over neighboring trees for possible sections of the path of rays
Simulate canopy PS, transpiration, light interception, etc. at every simulation interval
Initialize accumulators for cumulative water use, etc.
Loop over all time intervals
Retrieve weather data that has been particularized to this interval
Estimate eair and Tair inside the canopy, from previous time interval as a first estimate, to set environment of all leaves
Call subroutine (Ecalc) that computes canopy PS, transpiration (E) from leaf fluxes at all locations in canopy
Calculate errors (inconsistencies) between new e_{air} and T_{air} in the canopy, computed from canopy PS and E and aerodynamic conductance, and the initially estimated e_{air} and T_{air} in the canopy
Do 3 or more iterations of the same, until the errors are small enough
Extrapolate e_{air} and T_{air} to better estimate, using results of previous 3 iterations  call the subroutine named iterate to do this; also call subroutines econtrol, Tcontrol to keep estimated shifts in e_{air}, T_{air} from becoming large and causing instability in the iterations.
At this point, program has final estimate of wholecanopy PS and E
Report results and update accumulators
End of this simulation; allow a change in major physiological parameters and a rerun from start
Subroutine Ecalc  to compute canopytotal PS and E for one time interval
Declare variables and arrays
Set up common blocks (all the copious information that it needs enters this way)
Nested loops over all locations in the canopy
If not done in a prior run: compute light penetration probabilities to all these locations  else, retrieve the stored values
Compute total environment of leaves at this location  energy fluxes, especially
Initialize accumulators for contributions by leaves at different light levels at this location
Compute fractions of leaf area at this location that have irradiances in each of the preset ranges we are considering
Calculate PS and E rates of shaded leaves  requires a call to another subroutine (BBsolve) that gets the simultaneous solution of BallBerry equation of stomatal control and enzymekinetic equation of photosynthesis
Recalculate with 10% higher conductance  this is used to compute control coefficients (not of direct interest for present purposes)
Calculate
sunlit leaf PS and E
Loop over irradiance ranges on these (=loop over different orientations)
Calculate PS, E for midpoint irradiance in each range
Redo at 10% higher conductance, again for computing control coefficients
Return to calling program (main)
Subroutine BBsolve  to get the simultaneous solution of BallBerry equation of stomatal control and enzymekinetic equation of photosynthesis. Method is to find the calue of leafinternal CO_{2} (C_{i}) that solves both equations
Declare variables and arrays
Set up common blocks (all the copious information that it needs enters this way)
Get initial guess for maximal g_{s}  call other subroutines: Tleaf, which solves for leaf temperature (see below), and CiACssolve, which solves a quartic equation for C_{i} (see below)
Use g_{s} = b_{BB} as initial minimal estimate
For current guesses at gsmin, gsmax, also compute the midpoint, gsmid
For each of these points, compute A, hs, Cs and test if these satisfy the BallBerry equation
Reset gsmin and gsmax to the two points out of the three between which the sign of the error changes
Keep halving the interval until it is as small as the preset tolerance
Note: there are sections in this code that reset gsmin and gsmax if the range is too narrow (that is, there is no change in sign of the error between the 2 points)
Return to calling program
Subroutine Tleaf  solve for leaf temperature that satisfies energy balance (no net rate of energy gain)
Declare variables
Set up common blocks (all its information except the value of g_{s} enters this way)
Set initial estimate that T_{leaf }= T_{air}
Recall all other energy inputs
Compute rate of transpirational cooling from value of g_{s}, T_{leaf}, e_{air}, etc.
Iterate T_{leaf} until net energy gain is zero; use the NewtonRaphson method
Return to calling program
Subroutine CiACssolve  solves a quartic equation for the value of leafinternal CO_{2} partial pressure (C_{i}) that satisfies the enzymekinetic equation for photosynthetic rate and the transport equation for CO_{2}.
Declare variables
Set up common blocks (all information enters this way)
Calculate the enzymekinetic parameters of photosynthesis, which depend on temperature, using leaf temperature computed in subroutine Tleaf
Set up the terms in the quartic equation
Iterate for C_{i}, using the NewtonRaphson method  call the service subroutine F(C_{i}) to compute the error at each iteration
When the solution has converged, compute A and C_{s}.
Return to calling routine
Other subroutines:
F(Ci)  noted immediately above
julian  computes the Julian date for any specified month, day, and year
econtrol  limits the change allowed in canopy e_{air}_{ }on iterating for new value, to prevent numerical instabilities
Tcontrol  does same, for canopy T_{air}
iterate  uses the values of e_{air} and T_{air} from any 3 iterations, and their associated errors,
to compute a new estimate. It fits the
errors to a biquadratic in e_{air}
and T_{air}
Besides the source code, I have supplied compiled and linked versions that are directly executable. If you got these in an email, they will have phony extensions xxx, such as E7e.aragorn.xxx. I changed the extensions so that the files could be passed through email. Most email utilities filter out any attachments with .exe extensions, considering them rightly as potential malicious code. Simply change the name of the file so that the extension is exe. If you don't change it, you can still run the code, but you have to specify the full name, with the extension xxx.
The code is run in a command window. I keep one on my PC desktop as an icon. If you don't have this, then click on start à All programs à Accessories à Command prompt. In this window, you can start execution by simply typing the name of the file after the command prompt (which is usually the rightcaret, ">")  e.g., to run the complex model, type this:
> E7E.aragorn
You can run the code with manual input from the keyboard, entering each line of code. This is really impractical, when you have to specify the geometry of each of 81 trees. A practical way is to type all the inputs into a file and then use a redirect. Suppose you saved the input lines in the file pecan_input.txt. The command you type in is then:
> E7E.aragorn < pecan_input.txt
I have provided a complete sample file of the inputs, pecan_input.doc. It has (mostly) the values discussed above, and after each line there are comments that give the names of the inputs. You would have to save this as a simple text file, say, pecan_input.txt, to use it as input to the program.
The next consideration is that, by default, the output scrolls into the command window. There may be a huge amount of it. You can save all the output into a file for later perusal or later processing with another redirect. The command
> E7E.aragorn < pecan_input.txt > pecan_output.txt
will store the results in the file pecan_output.txt.
Interpreting and analyzing the output
There's some more work to be done here by me, to simplify and customize the output for our new purposes.
The essential elements in using the output are understanding the following items:
 what are the items and their units?
 what does one look for as patterns in the variables? à
 how does one distill the results to simpler quantities that we want, so that we can 1) design the experimental measurements the best way, to test the accuracy of our understanding, and 2) design the optimal strategy of deficit irrigation. For the 2^{nd} part, we want to extract some relatively simple quantities. The most critical quantity is, What is the shift in wateruse efficiency for a specified fractional reduction in transpiration? This needs to be computed for different stress levels (fractional reductions) and for different seasons (maybe just as VPD, which is another big controller of WUE). A description of the change in WUE with stress can be inserted directly into Ted's extant model of pecan growth and yield.
I note that WUE can be computed on many different bases. It is the ratio of carbon gain to water used. However, there are diverse measures of carbon gain, and diverse measures of water used. We have to be clear about which measures we are using, and we must use the final definition appropriate to our tests.
If we are doing leaf gasexchange measurements with a LICOR system, then carbon gain is the CO_{2} uptake rate in, usually, mol m^{2} s^{1}, and water used is leaf transpiration in the same units. A typical value might be 0.005.
If we want to compare our simulations with measurements done with eddy covariance, then we are comparing wholecanopy carbon exchange with wholecanopy evapotranspiration. Two more considerations enter. First, we have to account for all the respiration in the canopy and in the soil. This may be 5080% of leaf net photosynthesis during the day. We need a good model of trunk, root, and soil respiration to make this comparison. Then, we can relate net canopy CO_{2} flux, A_{can,net}, to the flux that is only in the leaves, A_{can,leaves}, and the total respiration per ground area, R_{can}:
_{}
This estimation is difficult. There is virtually a cottage industry in making these estimates. We can discuss this at length. Second, we need to account for the evaporation of water from the soil. This involves estimating the solar energy load on the soil, the boundarylayer conductance at the soil, and the soil wetness. This, too, involves sophisticated models, and I leave them out of my canopy models. This degrades their accuracy somewhat, because the fluxes of heat and water vapor from the soil alter the air temperature and humidity in the canopy.
 how do we find out which physiological and treestructural parameters are most important to measure accurately? That is, which ones show the greatest effect on nut yield and water use if they are varied within the bounds that are typical? The obvious way to quantify these sensitivities is to rerun the simulation with altered values of the parameters. For example, we may have entered 120 μmol m^{2} s^{1} as the photosynthetic capacity of leaves. We may have an estimate of the uncertainty as ±20 in the same units. We can then rerun the simulations at 100 and 140 as the value of PS capacity and check the final results, particularly the yield itself, and WUE and its stress response.
The quantities provided in the output of the complex model, and their units.
The model writes some data only to files, not the screen  particularly, the array of lightpenetration probabilities to all locations in the central tree, for every simulation interval.
Some diagnostics are printed before starting the simulation:
Light levels (PAR flux densities, in μmol m^{2} s^{1}) computed from total solar energy and inferred cloudiness, for 1^{st} 24 hours of weather data).
Keywords, prompts for changing parameter values (several places)
Diagnostics, for a few selected intervals (not of interest for current work; used for debugging in detail):
For interval # 121: Number of the iteration in which the iterations for canopy e_{air} and T_{air} converged.
Remanent errors in e_{air} and T_{air}_{ }after acceptable degree of convergence, scaled to tolerances for each and combined as random errors (square root of sum of squares)
For any interval, if e_{air} and T_{air} have not converged after 10 iterations: final errors in e_{air}_{ }(Pa) and T_{air} (ºC); on another line, the initial estimate and final 4 estimates of e_{air} (Pa); the last 4 estimates for stand transpiration rate (mol m^{2} s^{1}); the initial and final 4 estimates of T_{air} (ºC); the final 4 estimates for stand sensible heat flux (W m^{2}); the last 4 errors in e_{air} (Pa); the last 4 errors in T_{air} (ºC)
For interval #121: final estimates of e_{air} (Pa), T_{air} (ºC)in canopy, wholestand rates of transpiration (mol m^{2} s^{1}) and sensible heat flux (W m^{2}), final errors in e_{air} (Pa) and T_{air} (ºC), and absolute values of these errors.
For an interval that the user can specify: for many locations in the canopy: C_{i}/C_{a} (unitless, a measure of intrinsic wateruse efficiency), actual WUE (mol_{CO2} mol_{water}^{1}), leaf temperature (ºC), irradiance on leaves (μmol m^{2} s^{1}), surface relative humidity (unitless), and weighting of this leaf area (fraction of total leaf area it represents.
For EVERY interval of the simulation:
On one line: the number of the interval in the whole sequence (integer); the instantaneous transpiration rate of the whole tree (L/h); the instantaneous photosynthetic rate of the whole tree (mol_{CO2} s^{1}); the instantaneous control coefficient for g_{s}_{ }over transpiration rate and photosynthetic rate of the whole tree (fractional change in these rates per fractional change in g_{s}; unitless); initial and final values of e_{air} in the canopy (Pa); windspeed above the canopy (m/s); average irradiance on all leaves in the canopy (μmol m^{2} s^{1}); final and initial values of T_{air} in the canopy (ºC; note reversal of order); final errors in e_{air} (Pa) and T_{air} (ºC); tolerances used for e_{air}_{ }and T_{air}; number of the iteration in which e_{air} and T_{air} converged (integer).
On another line: running total of irradiances at all locations for this simulation interval (mol m^{2} s^{1}); same, accumulated over all intervals to date
At the end of the whole simulation period:
On one line: total leaf area on the tree (m^{2}); ground area dedicated to a single tree (m^{2}); average irradiance on all leaves (μmol m^{2} s^{1}); rate of total light interception by central tree during this interval (mol h^{1}); same rate, for soil below the crown; fraction of light intercepted by the tree, vs. the ground, in this interval (dimensionless); spacing of trees in EW direction (m); spacing of trees NS (m)
On another line: total water use by the central tree (liters); total leaf photosynthesis by the tree (mol_{CO2}, before respiratory losses); grand wateruse efficiency as their ratio (mol_{CO2} per L water).
On another line: control coefficients of gs over photosynthesis and transpiration, for the complete interval (unitless).
The quantities provided in the output of the simpler model, and their units.
For each day of the simulation, the model reports the following estimates:
* Unstressed case: leaf temperature (average in canopy), leaf internal CO_{2} partial pressure (C_{i,} in Pa), updated photosynthetic capacity of average leaf (V_{c,max}^{25}, in μmol m^{2} s^{1}), and stomatal conductance (g_{s}, in mol m^{2} s^{1}).
Same, for stressed case
* Unstressed case: wholestand transpiration rate (mm/d), and photosynthesis rate (g dry matter gained per m2 per d)
Same, for stressed case
* Unstressed case: nut mass increment (kg/ha) and updated leaf N content (percent)
Same, for stressed case
At the end of the season:
* Unstressed case: total water use (mm) and nut yield (kg/ha)
Same, for stressed case
* Unstressed case: final leaf N content (percent)
Same, for stressed case
Getting best estimates of both the key
input parameters and the detailed stress responses
Physiological parameters. Both the realworld behavior and the simulation results are likely to be most sensitive to three parameters  the photosynthetic capacity of leaves (Vcmax25) and the two parameters of stomatal control as captured in the BallBerry model (mBB, and bBB).
In an earlier section, I described how each is measured using gasexchange.
Vcmax25 for any chosen leaf is computable from photosynthetic rate under the current conditions, plus leaf temperature and leafinternal CO_{2} content reported (computed internally) by the LICOR system. This is expected to show variations of several kinds:
1) At any one
date in any given leaf, it should be stable over the day. However, it will vary from leaftoleaf,
primarily because the leaves see different average light levels. The plant has internal physiological controls
that allocate nitrogen  and, therefore, PS capacity  preferentially to leaves
that see the most sun and can do the most PS.
Earlier, I noted how Ülo Niinemets in
2) Seasonal changes in Vcmax25 are expected. A particularly strong cause should be the decline of leaf N content as the nutfill progresses, drawing N out to use in making nut protein. Therefore, we should make gasexchange measurement at three times during the growth season, at least. These need not be made on the same leaves each time; we simply want a sampling of the wholecanopy average.
3) Water stress may cause changes in Vcmax25. It has been shown a number of times that daily changes in water potential do not change Vcmax25, but stress over longer time appears to have an effect. I analyzed David Johnson's data from 6 dates in 2003 and found significant changes. Oddly, stress seemed correlated with higher values of Vcmax25, not lower. We can get some quantification of the stress effect from the same seasonal measurements proposed in item (2) immediately above. We will need to know the midday water potential of the leaves, and possibly the predawn value. Both can be obtained by using the pressure bomb method, which may be familiar to students. I have documentation on how to do these measurements, in any case, and similar documentation is freely available in the literature. One can avert the need to get up extremely early by covering leaves with light and watervaporimpermeable material such as aluminum foil (best painted white to keep it cool). The leaf does not transpire and, with some reservation, it remains near root water potential. It is practical to coordinate measurements of gas exchange and water potential on each leaf as one samples multiple leaves in a day.
Stomatal control parameters mBB and bBB are likewise computable from gasexchange measurements.
1) To get the parameter values that are in force on a given day and given leaf, we must deliberately sample a range of environmental conditions. Some discussion among all of us will be useful for setting up the best protocol. There are two basic ways to attain a range of conditions. One is to find leaves at different locations on the same date. They will vary in irradiances they are receiving, and also in temperature, but much less so in humidity and CO_{2} level. Another way is to deliberately impose a range of conditions on the same leaf, kept in the LICOR chamber. The LICOR is capable of enforcing higher or lower temperatures, using a Peltier cooler/heater. It can also change the humidity on the leaf (downward; with a lessthanportable dewpoint generator, it can set any desired humidity. With an optional light source, it can also enforce any desired light level (irradiance). Finally, it can control the CO_{2} level at the leaf, using small CO_{2} cartridges. This makes the protocol more complicated but may well be worth it. I note that the stomata of tree leaves respond rather slowly to these changes. They may stabilize at a new stomatal conductance only after 15 minutes for each change. That means that running, say, 4 different conditions will require somewhat longer than one hour, per leaf. One might get 5 or 6 leaves sampled in one day. This exercise is, nonetheless, important. We should plan on at least one day at each of three seasons doing these series of measurements.
2) At the same time that we are getting these measurements, we will see leaves that are at different water potentials, ψ (of the leaf itself, or the root system). We expect that greater water stress on the individual leaf, as lower ψ (more negative), will cause a reduction in the BallBerry slope, mBB. A form that several of us have used is
_{}
Here, m_{BB}^{0} is the unstressed value and β is a constant. We need not do any more sampling of leaves than is needed to estimate the effect of water stress on Vcmax25, described above.
3) There may be seasonal changes in mBB and bBB. We will see these in following the other parts of the protocol outlined above. There is no need for any separate sets of measurements. In David Johnson's data, the effect was pronounced; mBB almost doubled between premonsoon conditions (up to July) and post or duringmonsoon conditions (August and September).
Important soil variables to relate to stress and physiological responses to stress. Trees, like essentially all vascular plants, respond to water stress by changing stomatal conductance (g_{s}) and, in the longer term, the photosynthetic capacity of leaves (V_{c,max}^{25}). It is not yet clear if these changes are linked to leaf water potential (ψ_{L}) or to root water potential (ψ_{r}), or a mixture. In any event, we will want to measure these water potentials and then correlate them to shifts in the stomatal control parameters and in V_{c,max}^{25}.
We can infer root water potential, ψr, from soil water potential, ψs, and some additional information. There is a drop of water potential from bulk soil to roots that depends upon the rate of transpiration of the whole tree, E, and the hydraulic resistance from bulk soil to root surfaces (R_{soil}):
_{}
The soil resistance is a function of root length density (RLD), fineroot diameter, the volume covered by roots of the tree (determined from soil depth and radius of the root mass), and, finally, the soil hydraulic conductivity (k_{h}). This last conductivity is tied to the soil type and to the soil water potential by welltested relations.
Thus, measurement of ψ_{s} , the soil type, and some root parameters will give us both ψ_{s} and ψ_{r} that we need. Both Ted and I have programs to do the calculations. I have provided the source code and executable for my own Fortran program, testTsoil.f and testRsoil.xxx.
The soil water potential can be measured several ways. A direct method is with thermocouple psychrometers. This is a wellverified technology. There is a bit of work in calibrating individual psychrometers and in keeping them functioning in the field, as well a bit of care in placing them (e.g., not too close to the surface, where daily temperature variations can cause problems).
An indirect way is to measure soil volumetric water content, θ, and then use the soil moisture release curve for the given soil type to compute ψ_{s}. This can be satisfactory in some soils, but not in soils with significant clay content. In such soils, the change in ψ_{s} with changes in water content are very rapid. Small changes in the estimate of θ , from normal measurement error, generate large errors in estimates of ψ_{s} and even greater errors in estimates of ψ_{r}. To demonstrate this, I have made calculations at different soil water potentials (the most relevant measure). At each value of ψs, I computed the values of θ, Rsoil, and final root water potential, ψr, for a pecan tree transpiring 0.011 kg per s (equivalent to 0.5 mm/h for 9m diameter trees in an orchard, a moderately high rate but not maximal). I used root parameters that Ted has estimated for local pecans: a mass of fine roots per tree of 45 kg, a fineroot radius of 0.29 mm, a rooting depth of 1.2 m, and a rootmass radius of 4.5 m. I computed the rootlength density from the standard drymass density of roots, 250 kg m^{3}. I set the soil type to silty clay loam.
I computed the quantities noted above at 4 different soil water potentials of interest, from 0.1 to 0.4 MPa (1 to 4 bar). For each case, I also computed ψ_{s} and ψ_{r} (and R_{soil}) at water contents that were 1% higher or lower. Water content measured by timedomain reflectometry (TDR) has a typical error of this magnitude. This allows the estimation of the likely error in inferred ψ_{s} and inferred ψ_{r} with these errors in measuring θ. Here are the results:
ψsoil 
θ 
Rsoil 
ψroot 
ψsoil
with error in θ 
ψroot
with error in θ 
(Mpa) 
(%) 
(Pa (kg s1)1 
(Mpa) 


0.1 
22.27 
1.74E+06 
0.119 



23.27 
1.09E+06 

0.086 
0.098 

21.27 
2.92E+06 

0.115 
0.149 
0.2 
18.38 
1.63E+07 
0.380 



19.38 
8.57E+06 

0.169 
0.258 

17.38 
3.38E+07 

0.250 
0.621 
0.3 
16.63 
6.12E+07 
0.973 



17.63 
2.76E+07 

0.235 
0.538 

15.63 
1.50E+08 

0.395 
2.050 
0.4 
15.59 
1.57E+08 
2.125 



16.59 
6.32E+07 

0.303 
0.998 

14.59 
4.56E+08 

0.554 
5.566 
The errors in ψ_{s} in the region where stress at the root (ψ_{r}) becomes significant, around 1 MPa, are modest….but the errors in inferred ψ_{r} are enormous. At a true water potential of 0.4 MPa, a 1% measuring error with a TDR gives an error between +0.1 and 0.15 MPa (+1 to 1.5 bar), while the error in ψ_{r} is between +1.1 and 3.4 MPa!. Note that I have inflated soil resistance by a factor of 5, to account for the clumping of roots that reduces their efficacy in taking up water. This is the value used by F. Tardieu. It is difficult to measure accurately. If it is only 2 rather than 5, then the gradients between soil and root are only 40% as large, and so are the errors. The worst error at a soil water potential of 0.3 MPa would still be 0.095 MPa in ψ_{s} itself and 0.43 MPa in ψ_{r}.
The conclusion is that we must use psychrometry in this region. Even this will be challenging for estimating ψ_{r}. We will need to supplement it with measurements of leaf water potential on nontranspiring leaves.
Some other parameters, less variable. Some of the physiological parameters of pecans have not been measured with any systematic approach. One is the curvature of the light response of photosynthesis, thetaPS. However, it has a fairly small effect on the results, and the uncertainty in its value is likely quite small, given the studies on other plants. Another parameter is dark respiration of leaves, RdovrA. The same argument can be advanced here. Efforts to measure these parameters more accurately will not pay back any significant improvement in the results.
Some parameters of canopy structure will need to be measured. Certainly, we need the foliage density, fd, fairly accurately. Onetime measurements with a canopy analyzer after full leaf development should be sufficient. Leaf nitrogen content is also well worth measuring. It is related to PS capacity. Although we do not estimate PS capacity from leaf N in our study, it will be an effective way to do so in a decisionsupport system used at other sites. Moreover, for the optimization model, we will want to verify that leaf N declines as we have assumed, and especially how the decline differs between wellwatered and stressed trees.
Observing dynamic stress responses. The canopy models assume that the principal response to stress is a reduction in stomatal conductance, and, specifically, a reduction in the BallBerry slope, mBB. Another response is a lessened depletion of leaf N by nutfill. Certainly, the protocols described earlier will reveal these responses "in the large," that is, as seasonal trends. They will also reveal the intermediate mechanism, at least of the stomatal response. The reduction in water supply acts on stomata only by changing the root and leaf water potentials. We will be able to relate the changes in mBB to the water potentials and, we hope, to extract the proposed parameter β in the earlier equation. If this succeeds, we will be able to put more diagnostics into the decisionsupport system, such as leaf water potential under stress.
Scaling up the measurements of stress responses to whole canopies. The most relevant observations of stress responses would be measurements of wholetree or wholecanopy photosynthesis and transpiration. Unfortunately, there are no easy ways to do this. For single trees, transpirational water use all occurs throught the trunk, so that one might consider using sapflow measurements to quantify wholetree transpiration. My group has done this, with some success. We might consider this. It gives hourly or finer time resolution, but it requires that we know the sapwood area quite accurately. Recent research also shows that we cannot assume uniform sapflow at all azimuthal positions on the trunk  we might need several systems per tree; only about 12 gauges are available (they are very expensive), so we cannot sample many trees.
On a coarser time scale of a week or so, we can use soil water balance to get cumulative transpiration.
There are no good ways to measure photosynthesis on single trees. Researchers have bagged whole trees in a heroic effort, but it is very costly and timeconsuming, and the environment of the tree is changed so radically that it cannot be considered a good measure of stress responses.
For larger areas with many trees, covering at least 150 x 150 meters, one can use eddy covariance to get the fluxes of CO_{2} and water vapor in the air. None of our treatment areas is anywhere near large enough.
The conclusion is that we will have to rely on measurements done at the scale of leaves, plus models to scale up to canopies.
A
brief sketch of objective (3): Use
the models to parametrize an initial optimization model
The simple model ("big leaf") has already been used to estimate changes in photosynthesis, wateruse efficiency, and initial leafarea development under stress. The results, at only one stress level and the normal irrigation level, were used to parametrize the final results (PS, WUE, canopy LAI) as a function of irrigation fractions, f_{1} and f_{2}, during leafout and during nutfill. There is much work left, though it may not take a great deal of time:
* Validate the simple model against the more complex model. The simple model has to estimate the sunlit leaf fraction and the mean irradiance on leaves. I estimated the form of these functions but will need to verify them from output of the more complete model. I also used a single representative time of day. I will generalize the model to use a full diurnal course.
* Run either or both models at a range of stress levels. From this, generate continuous functions that describe, for example, how WUE varies with stress at any stress level. The complex model has no stress function in it; it will have to be introduced.
One unfortunate constraint in the experiments is that we cannot apply water stress during leafout, if only because it would require twice as many trees and twice the overall experimental effort. This means that we cannot get the responses of final leaf area to deficit irrigation during leafout. The simple model indicates that this is a crucial part of the whole management system. We will have to continue relying on estimates of the effect.
* Once experimental data are in hand, upgrade the models and redo these two steps.
How the models need to be modified
 a quick outline
Modifying the complex model to do the stress calculations will be done much the same as in the simple model. An input will let the user specify the stress level, as the fraction of normal irrigation water for this size of canopy. The program will execute a search in the value of mBB that achieves the desired reduction. I will need to introduce a leaf Nbalance module into the program, to account for withdrawal of N for nutfill; the same kind of calculation as in the simple model will suffice. As a simplification, the weather data will be modified to present the same weather each day  really, just a single day, characterized by a given vaporpressure deficit. (I believe, we will also need to see if air temperature also matters.)
In real operation of an orchard, water is not metered to the roots at exactly the same rate each day. Irrigation raises water use and the subsequent depletion of soil water generates some stomatal closure and reduced water use. The grower and the researcher can see cycles of changing transpiration rate and photosynthesis, and the underlying changes in stomatal conductance. These cycles will be resolved in a 2^{nd} version of the model. In this version, water inputs will be given as irrigation dates and amounts, with some precipitation. Soil water content at all times will be simulated using these inputs and with transpiration (and soil evaporation) as the output. This will be a new module, as a subroutine in the code. Soil water content will determine the soil water potential and its hydraulic resistance, in yet another new module (based on the program I have already made, testRsoil.f, recast as a subroutine). Canopy transpiration and the water potentials of roots and leaves will be set interactively. The water potentials will set new values of BallBerry slope and thus of transpiration. This will be a simple modification of the existing subroutines.
Training in the use of the models and use of equipment
Ted has asked me to lead some discussions and to offer some specific training in the use of models and of the LICOR gasexchange system.
I am deeply familiar with the design and operation of the LICOR systems, models 6200 and 6400. Ted has asked me to give 2 lectures in the SOIL 620 class, on 22 and 26 March 2009. In the past, I gave one lecture on theory and the other on handson operation. This year, I will give the theory lecture, but the 2^{nd} lecture will be on analyzing data from the instrument in order to estimate PS capacity and BallBerry parameters. We will use the data of David Johnson from 2003 for this exercise. Handson operation will be simulated as an online exercise, using software from LICOR. I will also show real field operation, outside of class, at dates to be negotiated.
I have used the canopy analyzer to get foliage density and can show others how to use it. The instrument is not in my laboratory but we can borrow it.
We will need to develop the simplest efficient method of actinometry to measure average light intensities at various locations in the canopy. We can discuss this shortly in arranged meetings.
For those who have not done pressurebomb measurements of leaf water potential, I can give demonstrations and also documentation.
It is reasonable to expect that the descriptions given here of the models and their use will not suffice to let other researchers use the models effectively for their own purposes, or even to understand the implications of the models' use in developing and interpreting experiments. I can offer sessions in which we use the models to explore their potential value.
E7E.aragorn.f and E7E.aragorn.exe (or .xxx)  the source code and the executable for the complex model. In its present form, it does not impose water stress. Its outputs are rather copious and perhaps excessive for the present purposes, and it uses units that agronomists typically do not use.
pecan_input.doc  a sample file of input data for the complex model. It is annotated. It must be resaved as a .txt file to be used for running the model.
ETstressed_full2.f  the simpler, bigleaf model for photosynthesis and transpiration in a pecan orchard, and their responses to water stress caused by supplying only a fraction of full water use.
ETstress.in.txt  a sample input data file to run the above.
computing_PS_capacity.xls  a spreadsheet that lets the user compute photosynthetic capacity from measured values of photosynthetic rate, leaf temperature, and internal CO_{2} partial pressure.
BB.xls  a more complete spreadsheet that uses more LICOR data and computes the BallBerry index, for datafitting that yields the BallBerry parameters mBB and bBB.
testRsoil.f and testRsoil.xxx  the source code and the executable for calculating soil water content, soiltoroot hydraulic resistance, and root water potential from soil water content, soil type, and root structure.