Climate change and the complex dynamics of green spruce aphid–spruce plantation interactions

Aphids can have a significant impact on the growth and commercial yield of spruce plantations. Here we develop a mechanistic deterministic mathematical model for the dynamics of the green spruce aphid (Elatobium abietum Walker) growing on Sitka spruce (Picea sitchensis (Bong.) Carr.). These grow in a northern British climate in managed plantations, with planting, thinning and a 60-year rotation. Aphid infestation rarely kills the tree but can reduce growth by up to 55%. We used the Edinburgh Forest Model (efm) to simulate spruce tree growth. The aphid sub-model is described in detail in an appendix. The only environmental variable which impacts immediately on aphid dynamics is air temperature which varies diurnally and seasonally. The efm variables that are directly significant for the aphid are leaf area and phloem nitrogen and carbon. Aphid population predictions include dying out, annual, biennual and other complex patterns, including chaos. Predicted impacts on plantation yield of managed forests can be large and variable, as has been observed; they are also much affected by temperature, CO2 concentration and other climate variables. However increased CO2 concentration appears to ameliorate the severity of the effects of increasing temperatures coupled to worsening aphid infestations on plantation yield.


1
The impact of climatic change on forests, pest populations, and their interactions, has 2 been the focus of much work over the last several decades. The biology is complex, and 3 involves both direct and indirect effects of multiple climate variables. Elevated CO 2 , 4 changing temperature, and their interactions have direct impacts on plant growth and 5 quality (from the herbivores' perspectives [1]). Many invertebrate pest species' 6 population dynamics are highly sensitive to ambient temperature means and 7 variances [2]. While a great deal of empirical experimental work has been done on 8 species of economic consequence, their long-term responses to climatic change can only 9 really be assessed with models. In this paper, we construct a mechanistic ecosystem 10 simulation model of the interaction between Sitka spruce trees (Picea sitchenis) and 11 their herbivores the green spruce aphid (Elatobium abietinum).

22
The aphid feeds on phloem sap [5]. It partially defoliates but rarely kills its host; it 23 can depress annual growth by 10-50% [6]. The impact of such infestations is generally 24 summarized by its effect on observed yield of stem wood during and over a rotation. A 25 rotation length of 60 years does not allow studies which directly address the problem to 26 be easily executed. Therefore research remains mostly empirical and short term. Randle 27 and Ludlow ([7] p. 33) state that "The ideal model for defoliation studies remains to be 28 developed" and this seems to remain largely true (but, for studies of aphids in other 29 systems, see Table A.1 for summary). 30 Interest in this particular tree-aphid system's response to climatic change dates back 31 to at least the mid-1990s. Straw [8] summarized the assessment at that time as follows 32 (p. 134): 33 "Defoliation of Sitka and Norway spruce by the green spruce aphid 34 (Elatobium abietinum) is limited in the UK primarily by periods of cold 35 weather which reduce the number of aphids overwintering. If winters 36 become milder, as current models of climate change predict, then the aphid 37 is likely to become more abundant and years with severe defoliation more 38 frequent. In such circumstances the productivity of spruce will decline." 39 1.2 Models of aphids and climatic change 40 Probably due to their economic importance, aphids have been the focus of numerous 41 modelling studies. These studies fall into three main types: statistical, agent-based, and 42 mechanistic. 43 44 In the context of climatic change, statistical models of insect responses are largely 45 so-called species distribution models (SDMs). Popular tools include Maxent (a form of 46 presence only logistic regression modelling), Genetic Algorithm for Rule-set Production 47 (GARP; [9]) and several others. For a comparsion of these techniques, see [10]. While 48 these methods have been popular for modelling the potential impact of climatic change 49 on non-aphid insects, they have rarely been used for aphids. There are plenty of   (Fig 3). The monthly data [38] were supplemented by monthly data from Clino [39]. 141 Daily values were obtained by linear interpolation of monthly data (see chapter 7 and 142 Section 7.5.3 of Thornley [34]). Daily data were applied to give diurnally changing data 143 for air temperature (Fig 3A), relative humidity (Fig 3A), radiation (which is calculated 144 from sunshine hours; Fig 3B), and wind (not shown). Soil temperature and rainfall were 145 assumed constant over each day. The source program, efm.csl, gives details (see 146 Appendix C).  Soil temperature was assumed to be diurnally constant and equal to mean daily air 163 temperature; it varies from a minimum of 1.45°C on 16 January to a maximum of 13.5 164°C on 16 July. Diurnal air temperature variation is a maximum of 5°C on 16 May and a 165 minimum of 2.15°C on 16 December. 166 Relative humidity (RH) was assumed to be at its daily maximum at dawn and at its 167 daily minimum at 15 h. The RH daily maximum has a minimum of 0.75 on 16 May 168 when the RH daily minimum is 0.65 and the RH daily minimum has a maximum of 0.88 169 on 16 December when the RH daily maximum is 0.91.   for incremental increases in temperature from −2°C to +3°C. 350 vpm CO 2 is shown in black, and 700 vpm CO 2 is shown in red. Notice the qualitative differences in population dynamics that emerge as temperture and CO 2 change. Next to each time series is the respective spectral density (arbitrary units) calculated from the last 32 years of each time series; the frequency (x-axis) has been re-scaled to display in years. Go to figure.

Statistical models
We can compare the dynamics depicted in Fig 4 to those seen in Fig 5. Here we plot 214 the normalized trap captures of alate E. abietinum using data extracted from Day et 215 al. [42]. We can see that this time series too contains periodic structure verging on 216 chaos. The spectral decomposition has less structure than we see in our simulations.

217
While we show the spectral decomposition for the same length of time (32 years) our 218 'sampling' of that period is 73 times more dense. Nevertheless, Fig 5 shows one feature 219 that is also seen in our simulations, within year cycles (see e.g. , Fig 4 +1°C, either vpm 220 CO 2 concentration).   Fig 3). The CO 2 concentration is denoted by the shading (dark color: 350 vpm, light color: 700 vpm). The aphid condition is denoted by the symbol (•: without aphids, ■: with aphids). Y C is the yield class (m 3 ha −1 y −1 ), defined as the volume (m 3 ) of timber harvested at the end of a rotation per hectare averaged over the duration of a rotation of 60 years; C sys (kg C m −2 ) is the total C in the system at the end of the 60 y rotation. Results are shown for the no-aphid infestation situation (a steady state as applied in Fig 4) and for that when aphids are applied at time zero [Eq (64)]. Go to figure.
5.3 C-sequestration (C sys ) and plantation yield (Y C )

234
The carbon-sequestration (C sys ) and plantation yields (Y C ) follow from the effects of 235 aphids density dynamics, and temperature and CO 2 concentration changes. Fig 7 shows 236 the four combinations of aphid presence and CO 2 for each temperature increment. It is 237 clear from the figure that doubling CO 2 results in greater C-sequestration and greater 238 plantation yields at temperatures of −1°C to +4°C. The presence of aphids results in 239 lower C-sequestration and lower plantation yields. In general these two effects seem to 240 be approximately additive. The dynamics of the aphid populations is interesting due to its variety, ranging from 243 aphids dying out (Fig 4A, −2°C) to various degrees of chaos with annual, biennial and 244 triennial cycles (Fig 4). It is not surprising that chaos occurs rather robustly in the 245 aphid-plant context. There is much non-linearity prevalent; there are far more than the 246 minimum two state variables required for chaos; the equations are non-autonomous; and 247 the intrinsic time scales of the dynamics of the aphid sub-model and the tree sub-model 248 are very different and incommensurate [43,44]. This chaotic variety indicates that 249 focusing on management/control strategies in order to minimize adverse consequences 250 may be difficult, even with the assistance of a process-based model (for further 251 consideration of chaotic aphid population dynamics see [45][46][47][48][49] [25][26][27][28][29] is in (a) the treatment (mechanistic or otherwise) of the impacts of rising 257 CO 2 on plant growth and (b) the use of the model to study impacts of climatic change. 258 Both the present study and Newman et al. found that the mechanistic consequences of 259 rising temperature tends to benefit aphid population growth while rising CO 2 tends to 260 be detrimental due to the changes in the plant caused by the rising CO 2 (see also [26] (Fig 7) and on leaf area index (Fig 6) and aphid dynamics (Fig 4). None of the previous 276 modelling exercises have been able to provide such an integrated picture of the possible 277 impacts of climatic change, let alone for plantation ecosystems.  for models which can be represented by ordinary differential equations using the 317 'rate-state' formalism (rate of change is a function of state variables + driving variables 318 + parameters). A significant advantage of ACSL is that it is non-procedural-that is, 319 the statements can be placed in any order; in the first step, translation, places the 320 statements in an executable order, which can then be compiled and linked as usual. The 321 non-procedural capability allows the programmer to put the statements in an order that 322 makes biological sense. unmanaged forests. The model is based on simplified physiology and biochemistry with 341 soil and water sub-models. The efm couples carbon (C), nitrogen (N) and water, fluxes 342 and pools and provides stoichiometric balancing of the items represented.

343
Here we used the efm in the evergreen forest mode, parameterized for Sitka spruce 344 growing in a north British environment (Section 4.1).

381
The phloem C:N substrate ratio [kg C (kg N) The aphid submodel is shown schematically in Fig 2. The model has 10 state variables. 434 There are two morphs of the aphid: alate (winged; ala) and aperterous (wingless, apt). 435 Each morph has four juvenile stages called instars, and an adult stage. All aphid state 436 variables have units of numbers of aphids per stem. All notation is summarized in 437 Appendix E. 438 We define the total number of alate instars (a lai ), alate adults (a laa ), the total alates 439 (a la ), the same for apterous aphids (a pti , a pta , a pt ) and the total number of aphids (a ph ) 440 as: At t = 0, a ph = 10 aphids per stem.

442
The spruce aphid does not colonize and feed on the current year's foliage [90]. From 443 the four foliage categories in the efm, only the last three are used here for the calculation 444 of aphid densities on foliage (A leaf,aph ). Dividing by the colonized leaf area per stem,

450
Using the aphid state variables, the C and N contents of the aphids per stem are 451 (units: kg aphid C, N stem −1 ): The C and N contents of each morphological aphid form (Fig 2)   There is a single input, from the developmental output flux of the fourth apterous instar 462 (a pt4 , Fig 2), O apt4→a , calculated in Eq (42) (the output is assumed to become an input 463 without loss). The input, I apta , requires C and N fluxes, of I Capt4→a and I Napt4→a : Units of the first of Eq (8)  the conclusions regarding the impacts of climatic change [29]. 484 It is assumed that aphid mortality, k aph,mort , depends on air temperature, T air and 485 phloem N, N phloem ; the specific rates (d −1 ) are calculated independently [Eqs (9) and 486 (10)] and then combined in Eq (11). Adult apterous mortality rate is then O apta→mort 487 (aphids stem −1 d −1 ), given by Eq (12).

488
The specific air-temperature-dependent mortality rate, k aph,Tmort (d −1 ), is given by 489 a skewed inverted parabola (Fig D.1): If air temperature T air equals the optimum temperature for minimum aphid mortality, 491 T aph,mort,opt , then mortality rate k aph,Tmort = k aph,Tmort,opt is zero. If temperature

492
T air = T aph,mort,0 (the second reference point, taken here to be 0°C), then The most significant part of the response drawn in Fig D.1 is the increase in 499 mortality as the temperature is lowered. This causes the aphid infestation to become 500 less severe or even zero as the temperature is lowered (Fig 4). It is assumed that the overall mortality rate (d −1 ), k aph,mort , is given by summing 512 the temperature-dependent and the nutritionally-dependent mortality rates [Eqs (9) and 513 (10); assuming aphid temperature = air temperature, T air (°C)]. Therefore: 514 k aph,mort = k aph,Tmort + k aph,Nmort The mortality rate varies between 2% and 6% per day, depending on temperature and N 515 nutrition (Fig D.3B).
With Eq (11) for specific aphid mortality rate, k aph,mort (d −1 ), output from the 527 apterous adult pool, a pta , to mortality is (aphids stem There are output (O) C and N fluxes associated with this aphid flux (using an obvious 529 notation): Units are kg C, N stem −1 d −1 . Aphid C and N contents are given in Eq (7). These 531 fluxes are input to the soil surface litter pools [Eqs (108) and (118)].

533
Foliage pruning, if applied, occurs at a rate of k le→prune (d −1 ). This gives rise to outputs 534 (losses) of aphids and associated fluxes of C and N (kg aphid C, N stem The foliage pruning constant, k le→prune , is not used in the simulations presented here 536 and it is set to zero. This aphid output flux is added to the mortality flux in Eq (15).

537
The total output of adult apterous aphids, O apta , is obtained by adding the The differential equation (apterous aphids stem −1 d −1 ) and initial value for the state 541 variable, a pta , for the apterous adult aphids is: a pta (t = 0) = 0 apterous adults stem −1 .
The input and output terms are given by Eqs (43) and (15). converted to wingless (apt) and winged (ala) 1st instar nymphs, a pt1 and a la1 (Fig 2). 547 These can have different N contents [Eq (7)]. Account must be taken of this before total 548 and fractional fecundities can be calculated [Section D.3.4; Eq (29)].
Here default values of the five parameters are given; the shape of the default response is 553 illustrated by the continuous lines in Fig D.4A-D. The default is a cubic. f (T ) is only 554 non zero between temperatures T 0 and T ′ 0 . Fig D.4A shows the effect of varying the 555 parameter q 1 . The response is initially sigmoidal if q 1 > 1. In Fig D.4B, q 2 is varied. It 556 can be seen that the steepness of the high-temperature decrease is highly dependent on 557 q 2 for q 2 < 1, although the early part of the curve is much less affected by q 2 . In Taking T = T air and with the default parameters given, Eq (17) is used to modify 564 many above-ground (shoot) plant processes (where it is reasonable to take shoot 565 temperature equal to air temperature, T air ), as well as the temperature dependence 566 assumed for some aphid processes, e.g. emigration [Eq (61)].

567
The effect of temperature on aphid fecundity, f Taph,fec , is obtained by using air 568 temperature, T air and the standard temperature function in Eq (17)  (21) f Taph,fec (T air ) is given in Eq (18). I Npapta,max20 is the maximum amount of N that can 584 be processed (p) at the reference temperature of 20°C. I Npapta,max is illustrated in Neither of the two arguments is affected by the presence of aphids-e.g. which will 591 depress N phloem (Fig D.5B). In fact the processing-limited maximum is about 592 eight-times as large as the volume-limited maximum (Fig D.5C This is illustrated in Fig D.5B and Fig D.5D.

598
The C intake per aphid (kg C aphid −1 d −1 ) and the N and C intakes per stem (kg N, 599 C stem −1 d −1 ) are respectively: The right side quantities are given in Eqs (23) The N ingested is assumed to be completely converted to wingless (apt) and winged 603 (ala) first instar nymphs, in pools a pt1 and a la1 (Fig 2). Because these two forms can This is drawn in Fig D.6. K Tapt is the half-maximal response temperature. q Tapt   617 determines the steepness of the response; q Tapt = 1 gives the familiar Michaelis-Menten 618 response. 619 D.3.2 Alate:apterous ratio as a function of aphid total density 620 This is ρ aph [aphids (m 2 leaf area) −1 ], Eq (5). Higher aphid densities cause more of the 621 apterous offspring to be alates (winged). As in Eq (25) and Fig D.6 above for relative 622 air temperature, a positive sigmoidal dependence on aphid density, ρ aph , is assumed: (ρ aph ) q ρapt + K ρaph,apt , q ρapt = 2, K ρaph,apt = 1000 aphids (m 2 leaf area) −1 .
(26) The half-maximum density point is where ρ aph = K ρaph,apt ; the value given is equivalent 624 to aphids being on a square grid 3 cm apart and is of the same order as the annual peak 625 in ρ aph .
The half-maximal phloem N concentration of K Napt = 4 kg N m −3 is equivalent to a 634 foliage N substrate concentration of N le = 0.04 kg N substrate (kg structural dry 635 matter) −1 [Eq (2)].

636
Combining these three factors multiplicatively [Eqs (25), (26) and (27)], therefore 637 the fraction of offspring from apterous adults which are alate (winged) and are apterous 638 are given by: These three components and the outcome for f apta→ala1 are illustrated in Fig D.7.

640
The relative contributions to the fraction of apterous offspring assigned to the alate 641 pathway are shown in Fig D.7 [Fig 2, Eq (28 (28). This is for an Eskdalemuir environment (Section 4.1) and the first year of growth of typical spruce plantation (Section 4.2) infected with ten alate adults at time zero [Eq (64)]. A, air temperature, T air , via Eq (25). B, aphid density, ρ aph , via Eq (26). C, phloem N, N phloem , via Eq (27). Here, the introduction of aphids at time zero, depresses N phloem , the lower continuous line, resulting in a higher fraction destined for alates (upper dashed line). D, the three factors shown in A, B, C are combined in Eq (28). Go to figure.

D.3.4 Total fecundity and the fractional fecundity of apterous adults 651
Now that we have calculated the fractions of offspring from apterous adults which are 652 alate and which are apterous, f apta→ala1 and f apta→apt1 in Eq (28), we are able to 653 calculate the total fecundity of apterous adults. It is assumed that N ingested by 654 apterous adult females per stem per day, I Nphloem→apta [Eq (24)], is converted without 655 loss into first-instar nymphs, giving a total rate of output to first instar offspring of The aphid N contents are given in Eq (7). Equation (29b) defines the total fractional 658 fecundity of apterous adults, f fecapta (d −1 ); this is the total output of nymphs from 659 apterous adults (a pta ) towards the alate and apterous first-instar pools divided by the 660 number of apterous adults (a pta ). The transfer of nymph outputs into inputs into the 661 first-instar pools is assumed to occur without loss giving (aphids stem

D.3.5 Associated N and C fluxes 663
There are output N fluxes from the a pta pool associated with these aphid fluxes of (kg 664 N stem The right-side quantities for the first two equations are in Eq (30) and Eq (7).

666
There are carbon fluxes from apterous adults to respiration, offspring and honeydew 667 [input to the surface litter metabolic pool, Eq (107)]. The carbon intake of the apterous 668 adults is I Cphloem→apta [Eq (24), kg C stem −1 d −1 ]. It is assumed that a fraction, 669 f Caph,resp , of this is respired, giving a respiratory flux from apterous adults of (kg C 670 stem −1 d −1 )
This value of 0.5 is typical of respiration losses [87].

674
The rest is honeydew (kg C stem (24), (32) and (33)] This output is input to the surface litter metabolic pool [Eq (107)].  The nutritional effect is small: a 10% decrease in development rate if N phloem [Eq

685
(2)] is zero so that f Naph,dev,min = 0.9. The nutritional modifier of developmental rate, 686 f Naph,dev , is assumed to be given by a Michaelis-Menten term with: f Naph,dev,min = 0.9, K Naph,dev = 1 kg N m 3 . (36b) The K value, giving half-maximal nutritional effect, is equivalent to foliage N substrate 688 concentration of N le = 0.01 kg N substrate (kg structural dry matter) Summing, over all τ apt1→a,20 = 15 d.

698
The overall value for the transition from 1st apterous instar to apterous adult agrees 699 reasonably with Dixon ([4], p. 109, table 6.1). Our attempts to relate present 700 assumptions to the data, equations and figures of Duffy et al. [63] were unsuccessful. It 701 is assumed that alates develop 15% more slowly than apterous aphids [Section D.9, Eqs 702 (9) and (80) As mentioned in the first paragraph of Section D.1.2, survival is a concept used by some 705 authors (e.g. [63]). In this subsection we show how it is related to the present 706 formulation.
The latter is drawn in Fig D.8B and the quantities determining this in Fig D.8A. It is convenient to deal with these pools together (Fig 2). The four instars have state 719 variables, a pt1,...,4 (aphids per stem). First the pool inputs are defined, then the outputs 720 and last the differential equations.
The inputs to the first pool are from the reproduction of both apterous and alate adults 724 [Fig 2; Eqs (30) and (74) The developmental rate constants (d −1 ) are given in Eq (38). The state variables, a pt , 733  = 1, . . . , 4 (aphids stem −1 ) are defined by Eq (47) (Fig 2). The outputs are transferred 734 without loss, providing the inputs (aphids stem −1 d −1 ) of I apt2−4 in Eq (41) and also 735 the input to apterous adults of Mortality. It is assumed that the same mortality rate constant, k aph,mort [Eq (11)], 737 applies to all four instars, leading to outputs of (aphids stem In the last two lines we calculate total instar mortality, O apti→mort and total apterous 739 aphid mortality, O apt→mort .

743
The last two of Eq (45) give the total apterous instar output to pruning and with the

747
Combining Eqs (42), (44) and (45), the total outputs from the four pools are (aphids 748 It is assumed that the N taken up from the phloem is the amount required to increase 755 the N component of mass of the aphids which enter the next pool. Therefore, using an 756 obvious notation (kg N stem −1 d −1 ): Aphid N contents are given in Eq (7); developmental outputs in Eq (42). Total N intake 758 from the phloem by the four apterous instars (apti) is (kg N stem −1 d −1 ): The C flux accompanying this N flux is (kg C stem −1 d −1 ): The phloem C:N ratio is given in Eq (3).

761
The C taken from the phloem is used for respiration, increase in instar mass and 762 honeydew. The earlier treatment is followed [Eqs (32)- (34)].

763
The amount respired is (kg C stem where fraction, f Caph,resp , is defined by Eq (32).

765
The increase in instar mass from aphid development (dev) also requires C from the 766 phloem. The four input (I) components of this and the total phloem C output to 767 supply these four inputs are (kg C stem Equations (42) and (7) are used in the above.

769
Excess C is excreted as honeydew (hon), which is input to the surface litter Eqs (50), (51) and (52) are used here.

775
Pruning similarly gives rise to C and N outputs of (kg C, N stem Equations (45) and (7) The total input to the alate adult pool is (aphids stem We have used Eq (7).  Mortality. The aphid flux to mortality, O alaa→mort (aphid number stem −1 d −1 ), is 791 assumed to depend on temperature and nutrition. The specific rate constant, k aph,mort 792 (d −1 ) assumed to be the same for all aphid morphs, is calculated in Eq (11). The output 793 from the adult alate pool to mortality and the accompanying C and N fluxes are [cf.
It is assumed that ten alate (winged) adult aphids fly in and colonize the stem at a The effect of temperature on alate fecundity is assumed to be the same as in Eq (18)   . The state variable for alate adults a alaa (Fig 2) is shown by the dashed-dot-dot-dot line [Eq (64)] with reference to the right side ordinate; initial value is 10 and is therefore not shown. Go to figure.
Here the total 1st instar progeny output from alate adults is O alaa→axx1 (aphids stem The N fluxes associated with these fluxes of aphids are (kg N stem The carbon intake of the alate adults is I Cphloem→alaa [Eq (69), (kg C stem There are carbon fluxes (kg C stem −1 d −1 ) to respiration, offspring and the residual C 866 flux goes to honeydew [cf. Eqs (32) to (34)]:

867
R alaa = f Caph,resp I Cphloem→alaa . (76) Here on the right side of Eqs (76) to (78) we have used Eqs (32), (69), (74)  There are four instars, with state variables a la1,...,4 (aphids per stem) (Fig 2). The 879 treatment is the same as above for the apterous instars (Section D.5) where the pools 880 are treated together, defining first the inputs, then the outputs and last the differential 881 equations. The equations are summarized below.
Rate constants are in Eq (79). State variables are defined in Eq (87). These outputs are 895 assumed to occur without loss, providing the inputs (aphids stem Eq (81) and also a developmental input to alate adults [see Eq (57)] of Mortality. These are [with Eqs (11) and (87); cf. Eq (44)] (aphids stem −1 d −1 ): It is assumed here that the same mortality rate constant applies to alate and apterous 899 instars. The last two of Eq (84) give the total alate instar and the total alate output to 900 mortality.

901
Pruning. These are [with Eqs (14) and (87); cf. Eq (45)] (aphids stem −1 d −1 ): The last three of Eq (85) give first, the total alate instar output to pruning, O alai→prune ; 903 next, with the first of Eq (60), total output of alate aphids to pruning, O ala→prune ; and 904 last, the total aphid output to pruning, O aph→prune , using the last of Eq (45 The C, N mortality fluxes above are captured in Eqs (108) These fluxes are captured in Eqs (110) and (120) and thence in Eqs (112)  Here, we calculate total mortality, emigration and fecundity for the aphids. Aphid fluxes to mortality are summed to give a total of O aph→mort (aphid number Here we have used Eqs (12), (44), (59) (59)]. In Fig D.9 the fractional output fluxes to mortality and to emigration are shown. 949 In the absence of pruning [Eq (63)] the two fractional variables defined in Eq (99) add 950 up to unity. Emigration peaks in late summer, when the adult alate population is 951 maximum (Fig D.9) and air temperature is also favourable (Fig 3A).

952
Because, as far as the local population is concerned, emigration is equivalent to 953 mortality, we also calculate the flux to 'mortality + emigration', which includes the 954 emigration flux from alate adults [O alaa→emi , Eq (62)]. Using an obvious notation, The terms Eq (98a) have units of aphids stem −1 d −1 ; in Eq (98b) we give the yearly 957 accumulations (aphids stem −1 y −1 ). In fact, the emigration flux is mostly small 958 compared to the total aphid mortality flux (Fig D.9). In Eq (98) and in the simulations 959 shown in Fig D.
The fractional fate of alate adults is drawn in Fig D. Total offspring production rate from apterous and alate adults (1st instar aphids stem Here we utilized Eqs (30)  The ten state variable differential-equation components of a ph (Fig 2) of Eq 4, namely 976 Eqs (64), (87), (16) and (47) are added to give da ph dt : Putting immigration, I imm→alaa = 0 [Eq (56)] and noting that alate adults (a laa ) are the 978 only morph able to emigrate, therefore Units are aphids stem  This has two purposes. First, it is a useful check that the model has been consistently 985 formulated mathematically and is programmed correctly. Second, in long-term 986 ecosystem studies, quite small fluxes are important, as they can lead to degradation, or 987 otherwise, of the ecosystem. These are now made explicit. There are two C inputs to the aphid sub-model, from the phloem and from possible 990 immigration of alates. The C input from the phloem is [Eqs (90), (50), (69) and (24) Units are kg C stem Units are kg C stem −1 d −1 (Eq (107a)) and kg C (m 2 ground) −1 d −1 (Eq (107b)).

1007
Mortality of aphids gives rise to C fluxes of [Eqs (59), (94), (12) and (54)] Units are as in Eq (107). This flux is divided equally between the surface litter cellulose 1009 and lignin pools [Eq (126) and Eq (127); figure 5.1 in [34]].   Fig 4. Aphid density dynamics. Shown are the time series of aphid densities (ρ aph ) for incremental increases in temperature from −2°C to 1113 +3°C. 350 vpm CO 2 is shown in black, and 700 vpm CO 2 is shown in red. Notice the qualitative differences in population dynamics that emerge 1114 as temperture and CO 2 change. Next to each time series is the respective spectral density (arbitrary units) calculated from the last 32 years of 1115 each time series; the frequency (x-axis) has been re-scaled to display in years. Return to text.  1, Fig 3). The CO 2 1125 concentration is denoted by the shading (dark color: 350 vpm, light color: 700 vpm).

1126
The aphid condition is denoted by the symbol (•: without aphids, ■: with aphids). Y C 1127 is the yield class (m 3 ha −1 y −1 ), defined as the volume (m 3 ) of timber harvested at the 1128 end of a rotation per hectare averaged over the duration of a rotation of 60 years; C sys 1129 (kg C m −2 ) is the total C in the system at the end of the 60 y rotation. Results are 1130 shown for the no-aphid infestation situation (a steady state as applied in Fig 4) and for 1131 that when aphids are applied at time zero [Eq (64)]. Return to text.