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 575-571-2269




Dates of preparation: Jan./Feb. 2009


Prepared on consulting contract to Prof. T. W. "Ted" Sammis, Plant and Environmental Sciences, New Mexico State University



Table of contents:


     Purpose of this document


     Overall goal

     A prologue on what we expect, and how and why we are combining models and experiments

          Expected tree responses

          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

          Critical elements of a model

          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 complex model

               The importance of canopy structure.

               Photosynthesis by the whole canopy

               The complete environment of the leaf determines photosynthesis

               Solving all these coupled equations

               Accounting for respiration

               The canopy affects its own environment

               Some processes are still left out

          The simple model

               A big leaf

               Solving the same equations

               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

               Canopy structure

               Controls over desired accuracy and over reuse of information

          Inputs to the simpler, big-leaf model

          Basic structure of the code for the complex model

          Running the code(s)

          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

               Physiological parameters

               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

          The tasks for the models

          How the models need to be modified - a quick outline

    Training in the use of the models and use of equipment

          Gas-exchange system

          Other field measurements

          Canopy models

     List of other files supplied



Purpose of this document


            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 decision-support 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 California, Texas, and New Mexico, and several other states that are not involved in the current project.  Production depends upon irrigation, but water supplies for irrigation  in the near term appears likely to be cut severely in California (15-50% of normal).  In the long term, both climate change and relentless population growth and attendant diversion of water to municipal and industrial growth will reduce irrigation water supplies permanently.

            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 season-total 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 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.


Overall goal:


             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, gs.  This cuts leaf transpiration almost in proportion - not quite as much, because leaf cooling is reduced, and the rise in temperature raises the leaf-to-air gradient in water-vapor pressure.  The reduction in gs 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 CO2.  Consequently, water-use 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 long-lasting 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 CO2 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 - Berry model (1980 ff.) puts all this into a surprisingly simple mathematical form.  Another robust pattern is in stomatal control, the physiological setting of gs by light level, air temperature, CO2, air pressure, humidity, windspeed, and water stress.  Researchers developed huge numbers of expressions for the roles of each factor, but they were largely unsatisfactory outside of a limited range of conditions until Tim Ball and colleagues developed the Ball - Berry equation (Ball et al., 1987).  It is also deceptively simple (see later in this document), while incorporating a welter of verifed physiological feedback and feedforward controls.

            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 real-time responses by growers, that is, decision-support 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 gs 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 Ball-Berry equation.  If this truly captures most of the variance in gs 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 Ball-Berry 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 Ball-Berry 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 enzyme-kinetic equation of photosynthesis (the Farquhar - von Caemmerer - Berry equation), the equation for stomatal control (Ball - Berry), the transport equation for transpiration, etc.  We can never account for every single process - e.g., detailed biochemical responses to transient changes in lighting as leaves flutter.  We have to be judicious in selecting the most important processes.  We can also use simple, perhaps fixed descriptions of other processes that are "a bit" important - e.g., soil evaporation under a closed canopy, described as a fixed rate, an average over wet and dry times.  This is one reason that modelling is an art, though an art that always bears the test of science, not art.  It matters not that a model is elegant, if it is wrong.


            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 non-transferrability 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 "data-hungry," 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 root-length 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 poorly-known 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 "big-leaf" model described later) was constructed so that I could force the transpiration rate to be one-half the rate of a fully-watered orchard.  In particular, this was achieved by finding, by simple search, the magnitude of the Ball-Berry slope parameter, mBB,  that would do this.  This choice of mBB  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 n1 and n2, 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 whole-canopy transpirational rate (to a fraction f1 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 f2 of normal) decreases photosynthesis per unit leaf area, as the leading effect.

            * Water use rates per unit leaf area (E1 and E2) differ between the two intervals.   Weather during leafout is hotter and drier, so that E1 is higher than E2.


            From these premises, we can compute estimates of:

        * Total water use, as a function of the chosen irrigation fractions f1 and f2.

        * Total photosynthesis, hence, total nut yield, as a function of f1 and f2.


            Finally, we can find the optimal strategy of deficit irrigation:

            * We can normalize these by the total water use with full irrigation.  Call this FW, 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 FY.

            * Let's set FW to a desired (or, really,  required) value, say, ½ .

            * Then, for any choice of the fractional irrigation rate in interval 1, f1, there is a matching value of f2 that makes FW = 0.5.

            * For these values of f1 and f2, we can compute (from the earlier premises) the value of FY.

            * Finally, we can find the best value of f1 (and, hence, also f2), as the one that makes relative yield, FY, 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 f1 in, say, a spreadsheet (which I provide, shortly), or a Fortran program, and seeing for which value we get the highest FY. 


     We will now make an estimate, putting some quantitative expressions into the premises:

       * Leafout lasts two months (n1 = 2), nutfill lasts 4 months (n2 = 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 n1 that is lower, say, n1,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 n2, similarly.  I propose n2,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 less-than-proportional drop in photosynthesis).  Let's make a linear estimate, that the fractional cover (fractional leaf area development), fcover, responds linearly to the irrigation fraction, f1:


       We can easily set a and b from fcover reaching 1 when f1 = 1 (so, a+b=1) and from fcover reaching 0.7 when f1 = 0.5 (that is, 0.7 = a*0.5+b = a*0.5+[1-a]).  The last equation gives us a = 0.6 and b = 0.4.

   * Increase in water-use 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 WUE0, we get


       This form gives us the proper ratio at f2 = 1 and f2 = 0.

    * To get the photosynthetic rate (per leaf area) at reduced irrigation, we multiply transpiration rate (also per leaf area) by water-use efficiency.  We have defined the stress factor f2 as the relative rate of transpiration per leaf area.  This gives us the photosynthetic rate per unit leaf area, A, as


      * To get the whole-canopy 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*f1+0.4) from the earlier note.  The complete factor of proportion is then


       * This captures the effect of deficit-irrigation fractions, f1 and f2, 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 WR1), multiplied by the relative reduction in canopy transpiration, multiplied by the effective duration of leafout.  The relative reduction in whole-canopy transpiration is the relative reduction in rate per leaf area (simply f1) 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 WR1 is 1.5 times larger than WR2, we can write the total water use, W, as


       The rate when full irrigation is applied is found by plugging in f1 = 1 = f2.  Call this W0:


       The relative reduction in water is the total use with deficit irrigation, divided by the total use with full irrigation:


       The factors of WR1 cancelled out between numerator and denominator.

       *  As noted earlier, we are given the required cut in FW - say, 0.5.  Now we can choose f1 and see what f2 has to be.  For example, if f1=0.5, then f2 has to satisfy


       We can do the reverse, too - we can specify f2 and ask what f1 has to be to reduce total water use by half. 

       * The final task is finding, for any choice of f1 and thus of f2 (or vice versa), what is the effect on yield, as the value of FY…and then testing various values of the f's to get the highest FY…for a specified cut in total water use.  Here is a table from the spreadsheet  I made:


































































            In the last line, I made a rough optimization.  I saw that the best FY occurred between f2 = 0.9 and f2 = 0.8.  At the intermediate value of f2 = 0.85, we require f1 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 f2 = 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 fcover = 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 (f1 = 1), then we'd have to cut F2 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 decision-support 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 water-production function of pecan orchards and of how their water-use efficiency varies with stress level.  Several auxiliary tasks include generating descriptions of N-stress 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.


      The complex model

       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 Pdir that the direct beam arrives at full intensity and a probability 1-Pdir that is completely blocked at this location.  The diffuse beam arrives determinstically, at a fraction Pdiff 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, Aleaf,  depends not only on the irradiance (in photosynthetically active radiation between 400 and 700 nm in wavelength), but also on temperature, humidity, CO2 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 Aleaf….and, at the same time, the leaf transpiration rate, Eleaf, and the stomatal conductance, gs:

      1) The Farquhar - von Caemmerer - Berry model of Aleaf in terms of basic photosynthetic capacity (Vc,max, related to content of Rubisco enzyme, in essence, and closely related to leaf N content), CO2 partial pressure at the Rubisco sites (Ci), and leaf temperature. 


Here, Γ and KCO are temperature-dependent 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 near-infrared 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 boundary-layer 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 energy-balance equation iteratively.  The boundary-layer 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 Ball-Berry equation for stomatal conductance.  This was also mentioned above.  It takes the simple form


Here, A is the leaf PS rate, Aleaf, and hs and Cs are the relative humidity and the CO2 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 hs and Cs depend on A and E of the leaf and on the stomatal and boundary-layer conductances, using readily-derived formulas of physical transport phenomena.

            4) The equations of transport for CO2 and water vapor: there must be a drop in CO2 partial pressure from outside to inside the leaf, to support the inflow of CO2 during photosynthesis.  There must be a rise in water vapor, to support the transpiration rate.  This gives us the relation for internal CO2 (Ci, the quantity need by the PS enzyme-kinetic equation) to external CO2 partial pressure, Ca:


Here, Pair is the total air pressure (this comes in because we use so-called molar units of conductance, the modern convention, for various good reasons), and gtot,CO2 is the total conductance for CO2 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 gs until all the equations are solved simultaneously.  Consider a guess at the value for gs. The energy-balance equation has all the other quantities specified, so this lets one calculate the leaf temperature (with quick iterations).  We combine the enzyme-kinetic equation (1), with its parameters corrected for the leaf temperature, and the transport equation (4) to get a single equation for Ci.  When we use the form of the enzyme-kinetic equation generalized to handle light-limited photosynthesis, this becomes a quartic equation.  We solve it rapidly by iteration.  Now we have both Ci and A.  Finally, we rewrite the Ball-Berry equation to highlight the error in the solution, as


When we have the right guess for gs, F becomes zero.  We home in on the proper value of gs by a binary search.  We guess the min and max values that gs could lie between.  We compute F at each end, and then for gs in the middle.  The solution has to lie between the values of gs 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 m-2 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 two-week-average air temperature, Tmean.   We input the latter and calculate the respiration rate for any leaf, applying an exponential factor in actual leaf temperature, exp(0.07*(T-Tmean)).

            The rate of photosynthesis is not to be compared with net CO2 exchange of an orchard, because respiratory losses of CO2 (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 CO2 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 CO2 level, and convective energy transfer alters the air temperature in the canopy.  We have to iterate the solution for, particularly, the air temperature, Tair, and water-vapor partial pressure, eair.  At each iteration, we get a new eair and a new Tair…and then new canopy totals of A and E…which gives us new eair and Tair.  The iterations are prone to oscillate and diverge, so that I apply a limiting factor.  The changes in eair and Tair, from their values in "free" air above the canopy, depend on the boundary-layer (or aerodynamic) conductance of the canopy as a whole.  This depends inversely on windspeed, with a constant of proportionality that depends on canopy leaf-area 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 whole-tree 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 whole-tree transpiration, gives us the drop in water potential at the roots.  The root water potential may be used to estimate the drop in Ball-Berry slope, m.  In a second version of the model, we simply enforce a halving (or other change) in whole-tree transpiration.  We search for the value of m that achieves this.  There are consequent changes in leaf temperature, gs, and A, all of which we use to see how whole-tree 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 decision-support 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.


    The simple model

            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 big-leaf 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 leaf-area 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 randomly-oriented 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, gs, and E per unit leaf area using the same set of 4 simultaneous equations as used in the complex model  - that is, the enzyme-kinetic 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, Vc,max25, 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 whole-tree photosynthesis to the rate of growth of nuts.  It assumes that a constant fraction of photosynthate (newly-made 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 (f1 or f2) noted in the description of optimal deficit irrigation, above.  It does so, again as noted earlier, by searching for the value of the Ball-Berry slope that forces transpiration to be half of the unstressed (fully-watered) 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 leaf-area 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 decision-support 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 tree-structural 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 user-friendly 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 edited-down entries from the standard data on the site, 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 2-digit integer - with a leading zero if it is less than 10

    day, as a 2-digit integer

    hour, as a 2-digit 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 two-digit 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 3-digit integer)

    insolation (solar energy falling on a horizontal surface), in MJ per m2 per h, as a real number with two leading digits and 2 digits after the decimal point

   * date0wea - the date on which the weather-data file begins, which is not necessarily the date on which the simulation begins (the user may have chosen a convenient whole-month block for the weather file, for example) - as text between quotes, such as "07-01-03".  Hyphens are required between the parts of the date, which must be in month-day-year order.  Each element must be a 2-digit integer - e.g., 07 for July, not simply 7.

   * Ca - the partial pressure of CO2 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 CO2 in air (today, about 385 ppm; it varies notably only during inversions).  Example: the air pressure in Las Cruces is typically 88,500 Pa (also not measured at most stations; it varies between highs and lows by several percent; we have no hurricanes here to make major excursions).  Then, Ca = 88500 Pa * 385x10-6 =  34.1 Pa

   * 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) CO2 fixed per m2 per s, the same units reported by all common measurement systems, such as the LI-COR 6400.  This is the rate that photosynthesis would attain at high (non-limiting) light intensity and high (non-limiting) CO2.  It is used in computing the PS rate at any light and CO2 levels, as seen in earlier equations.  It is in direct proportion to the amount of Rubisco enzyme per unit leaf area, Rubisco being the rate-limiting 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 College Ave.  Capacity can be computed from data in any standard gas-exchange system used on leaves - the measured photosynthetic rate (at high light levels, note - say, greater than 1000 μmol m-2 s-1), the internal CO2 partial pressure (Ci), and the leaf temperature.  One simply inverts the equation given by Farquhar, von Caemmerer, and Berry, earlier.  Note that one must convert the value of Ci reported by the LI-COR from ppm as a mixing ratio to partial pressure as Pa.  The LI-COR reports total air pressure as millibars, so the final equation to get Ci in Pa is


Leaf temperature, TL, in ºC, is used to compute the two enzyme-kinetic parameters, Γ and KCO, 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 light-levels in the long term - say, from different positions in the canopy.  Ülo Niinemets in Estonia has done an impressive series of studies that cover this, finding that Vcmax25 scales closely with average light level.  This can be incorporated into the model by some clever, simple means (which I haven't done yet.)

   * thetaPS - curvature in the transition between light-limited and light-saturated PS rates.  At low light levels, the PS rate, A, is proportional to the light intensity, IL (irradiance on the leaves, in the PAR region of the spectrum):


Here, the subscript "LL" denotes "light-limited." At high light levels, PS is light-saturated, at the rate expressed by the Farquhar - von Caemmerer - Berry rate, above.  The transition is not sharp, that A is at the light-limited rate until it reaches the saturated rate and then suddenly is constant.  A very good representation of the transition is


Here, Amax is the light-saturated 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 curve-fitting  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 gas-exchange system is the net rate, after respiration.  We formulate the net rate of PS, A = Anet as the gross rate minus respiration, Agross - Rd, and then as Agross*(1-Rd/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 long-term (2-week) 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 Ball-Berry relation for stomatal conductance, as a dimensionless number.  This is constant over the medium term, but it may vary seasonally; we need some gas-exchange measurements over the seasons to determine this.  A mean value is about 10.  The value has a significant effect on PS rate and water-use efficiency. 

            The value of mBB can be extracted from a series of gas-exchange 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 LI-COR 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 Ball-Berry index, A hs/Cs, from the measurements at each condition, and then plots the values against the values of stomatal conductance, gs, at each condition.  The value of A comes directly from the LI-COR.  The relative humidity at the leaf surface beneath the boundary layer, hs, is offset from the value in free air.  It involves the boundary-layer (gbL) and stomatal conductances and the water-vapor partial pressures in free air (eair) and in the leaf (eL, computed from leaf temperature).  The equation is


            Another equation relates Cs, the mixing ratio for CO2 in air at the leaf surface, to the mixing ratio in free air, Ca (remember, not in Pascals), the photosynthetic rate, and the boundary-layer conductance:


These calculations are not programmed into the LI-COR.  One must import the data from the LI-COR 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 LI-COR runs.  It also calculates Vcmax25.  It uses data that David Johnson acquired in his M.S. research in 2003.  The many columns of LI-COR 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 Ball-Berry relation, in mol m-2 s-1.  This is obtained in the same data-fitting 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 boundary-layer 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 m2 per m3, that is in m-1.  This can be measured by various instruments, including a LI-COR LAI-2000 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)πR3 = 382 m3, a total leaf area of 382*0.6 m2 = 229 m2, and a ground area per tree of 9*9 = 81 m2.  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 well-documented model called SiB, developed by Piers Sellers and others.  The formula is gb,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 (E-W and N-S) 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 positive-x 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 semi-axes 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 more-vertical axis is tilting off-vertical 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 canopy-total photosynthesis much.

   * I00start - the maximal intensity of the direct solar beam, in μmol m-2 s-1.  In the clear air of Las Cruces, this has a value near 2000.  It is affected by solar angle, but not as sharply as a cosine function (which is applied also to the total flux density).  It is affected by the air mass that the beam must pass through; the program computes this correction.  The program also corrects, in a first approximation, for haziness and cloudiness.  It computes the total solar flux density in the PAR, as a fraction (1/2) of the total solar energy, converted to photon flux density.  If this density is less than the diffuse density plus the direct-beam intensity multiplied by the cosine of the solar zenith angle, the program corrects both densities, using an algorithm that I developed.

   * 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 light-penetration 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 light-penetration 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 Ball-Berry equation in computing stomatal conductance, gs, including the factor of surface relative humidity, hs.  If it is positive, then hs will replaced by hs0 in the calculation of gs.  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 1024-fold 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 light-penetration 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 9-m 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 10-minute 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 month-day-year, e.g., 07-01-03 07-31-03. 

   * 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 per-tree values to the rates per unit ground area (say, kg water per m2 per hour or mm per day).  For nonuniformly spaced trees, these parameters are meaningless.


Inputs to the simpler, big-leaf 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 gs.


            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 (non-executable 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 eair and Tair in the canopy, computed from canopy PS and E and aerodynamic conductance, and the initially estimated eair and Tair in the canopy

      Do 3 or more iterations of the same, until the errors are small enough

          Extrapolate eair and Tair 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 eair, Tair from becoming large and causing instability in the iterations.

      At this point, program has final estimate of whole-canopy 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 canopy-total 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 Ball-Berry equation of stomatal control and enzyme-kinetic 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 Ball-Berry equation of stomatal control and enzyme-kinetic equation of photosynthesis.  Method is to find the calue of leaf-internal CO2 (Ci) 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 gs - call other subroutines: Tleaf, which solves for leaf temperature (see below), and CiACssolve, which solves a quartic equation for Ci (see below)

      Use gs = bBB 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 Ball-Berry 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 gs enters this way)

      Set initial estimate that Tleaf = Tair

      Recall all other energy inputs

      Compute rate of transpirational cooling from value of gs, Tleaf, eair, etc.

      Iterate Tleaf until net energy gain is zero; use the Newton-Raphson method

   Return to calling program


   Subroutine CiACssolve - solves a quartic equation for the value of leaf-internal CO2 partial pressure (Ci) that satisfies the enzyme-kinetic equation for photosynthetic rate and the transport equation for CO2.

      Declare variables

      Set up common blocks (all information enters this way)

      Calculate the enzyme-kinetic parameters of photosynthesis, which depend on temperature, using leaf temperature computed in subroutine Tleaf

      Set up the terms in the quartic equation

      Iterate for Ci, using the Newton-Raphson method - call the service subroutine F(Ci) to compute the error at each iteration

      When the solution has converged, compute A and Cs.

    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 eair on iterating for new value, to prevent numerical instabilities

   Tcontrol - does same, for canopy Tair

           iterate - uses the values of eair and Tair from any 3 iterations, and their associated errors, to compute a new estimate.  It fits the errors to a bi-quadratic in eair and Tair


Running the code(s)


            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 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 right-caret, ">") - 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 2nd part, we want to extract some relatively simple quantities.  The most critical quantity is, What is the shift in water-use 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 gas-exchange measurements with a LI-COR system, then carbon gain is the CO2 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 whole-canopy carbon exchange with whole-canopy evapotranspiration. Two more considerations enter.  First, we have to account for all the respiration in the canopy and in the soil.  This may be 50-80% 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 CO2 flux, Acan,net, to the flux that is only in the leaves, Acan,leaves, and the total respiration per ground area, Rcan:


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 boundary-layer 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 tree-structural 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 light-penetration 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 1st 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 eair and Tair converged.

     Remanent errors in eair and Tair 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 eair and Tair have not converged after 10 iterations: final errors in eair (Pa) and Tair (ºC); on another line, the initial estimate and final 4 estimates of eair (Pa); the last 4 estimates for stand transpiration rate (mol m-2 s-1); the initial and final 4 estimates of Tair (ºC); the final 4 estimates for stand sensible heat flux (W m-2); the last 4 errors in eair (Pa); the last 4 errors in Tair (ºC)

     For interval #121: final estimates of eair (Pa), Tair (ºC)in canopy, whole-stand rates of transpiration (mol m-2 s-1) and sensible heat flux (W m-2), final errors in eair (Pa) and Tair (ºC), and absolute values of these errors.

     For an interval that the user can specify: for many locations in the canopy: Ci/Ca (unitless, a measure of intrinsic water-use efficiency), actual WUE (molCO2 molwater-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 (molCO2 s-1); the instantaneous control coefficient for gs over transpiration rate and photosynthetic rate of the whole tree (fractional change in these rates per fractional change in gs; unitless); initial and final values of eair 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 Tair in the canopy (ºC; note reversal of order); final errors in eair (Pa) and Tair (ºC); tolerances used for eair and Tair; number of the iteration in which eair and Tair 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 (m2); ground area dedicated to a single tree (m2); 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 E-W direction (m); spacing of trees N-S (m)

    On another line: total water use by the central tree (liters); total leaf photosynthesis by the tree (molCO2, before respiratory losses); grand water-use efficiency as their ratio (molCO2 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 CO2 partial pressure (Ci, in Pa), updated photosynthetic capacity of average leaf (Vc,max25, in μmol m-2 s-1), and stomatal conductance (gs, in mol m-2 s-1).

      Same, for stressed case

   * Unstressed case: whole-stand 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


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


            Physiological parameters.  Both the real-world 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 Ball-Berry model (mBB, and bBB).

In an earlier section, I described how each is measured using gas-exchange.

            Vcmax25 for any chosen leaf is computable from photosynthetic rate under the current conditions, plus leaf temperature and leaf-internal CO2 content reported (computed internally) by the LI-COR 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 leaf-to-leaf, 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 Estonia has done an impressive series of studies that cover this, finding that Vcmax25 scales closely with average light level.   We should do the same on our trees, selecting leaves that sample a range of average light levels, which we might call Iavg.  We can't estimate Iavg on an instantaneous basis, obviously.  We need to install some sensor that records a longer-term average, say, over several days or a week.  Then, we can go back to the same leaves, which have been tagged, and do gas-exchange measurements on them.  What kind of sensors can do this?  An accurate system uses small photodiodes clipped onto leaves and tied to a datalogger.  I have developed such a system.  Installation can be quite tedious, however.  A simpler system is actinometry.  One places a small patch of a weakly photosensitive material on the leaf.  Over time, sunlight causes a cumulative photochemical reaction.  Development of the material reveals the value of Iavg, given a calibration set of materials in full sun.

     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 gas-exchange 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 whole-canopy 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 pre-dawn 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 water-vapor-impermeable 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 gas-exchange 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 CO2 level.  Another way is to deliberately impose a range of conditions on the same leaf, kept in the LI-COR chamber.  The LI-COR 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 less-than-portable dew-point 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 CO2 level at the leaf, using small CO2 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 Ball-Berry slope, mBB.  A form that several of us have used is


Here, mBB0 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 pre-monsoon conditions (up to July) and post- or during-monsoon 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 (gs) and, in the longer term, the photosynthetic capacity of leaves (Vc,max25). 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 Vc,max25. 

            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 (Rsoil):


The soil resistance is a function of root length density (RLD), fine-root 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 (kh).  This last conductivity is tied to the soil type and to the soil water potential by well-tested 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

            The soil water potential can be measured several ways.  A direct method is with thermocouple psychrometers.  This is a well-verified 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 9-m 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 fine-root radius of 0.29 mm, a rooting depth of 1.2 m, and a root-mass radius of 4.5 m.  I computed the root-length density from the standard dry-mass 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 Rsoil) at water contents that were 1% higher or lower.  Water content measured by time-domain 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 with error in θ

ψroot with error in θ



   (Pa (kg s-1)-1












































































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.  One-time 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 decision-support 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 well-watered 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 Ball-Berry 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 decision-support 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 whole-tree or whole-canopy 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 whole-tree 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 time-consuming, 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 CO2 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 tasks for the models


            The simple model ("big leaf") has already been used to estimate changes in photosynthesis, water-use efficiency, and initial leaf-area 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, f1 and f2, 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 N-balance 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 vapor-pressure 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 2nd 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 Ball-Berry 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 LI-COR gas-exchange system.


Gas-exchange system


            I am deeply familiar with the design and operation of the LI-COR 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 hands-on operation.  This year, I will give the theory lecture, but the 2nd lecture will be on analyzing data from the instrument in order to estimate PS capacity and Ball-Berry parameters.  We will use the data of David Johnson from 2003 for this exercise.  Hands-on operation will be simulated as an on-line exercise, using software from LI-COR.  I will also show real field operation, outside of class, at dates to be negotiated.


Other field measurements


            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 pressure-bomb measurements of leaf water potential, I can give demonstrations and also documentation.


Canopy models


            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.



List of other files supplied


   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, big-leaf 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. - 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 CO2 partial pressure.

   BB.xls - a more complete spreadsheet that uses more LI-COR data and computes the Ball-Berry index, for data-fitting that yields the Ball-Berry parameters mBB and bBB.

   testRsoil.f and - the source code and the executable for calculating soil water content, soil-to-root hydraulic resistance, and root water potential from soil water content, soil type, and root structure.