Expected Impacts of Mixing European Beech with Silver Fir on Regional Air Quality and Radiation Balance

Abstract

The anticipated climate change during the next decades is posing crucial challenges to ecosystems. In order to decrease the vulnerability of forests, introducing tree species’ mixtures are a viable strategy, with deep-rooting native Silver fir (Abies alba) being a primary candidate for admixture into current pure stands of European beech (Fagus sylvatica) especially in mountainous areas. Such a change in forest structure also has effects on the regional scale, which, however, have been seldomly quantified. Therefore, we measured and modeled radiative balance and air chemistry impacts of admixing Silver fir to European beech stands, including changes in biogenic volatile organic compound emissions. An increased fraction of Silver fir caused a smaller albedo and a (simulated) larger evapotranspiration, leading to a dryer and warmer forest. While isoprene emission was negligible for both species, sesquiterpene and monoterpene emissions were larger for fir than for beech. From these differences, we derived that ozone concentration as well as secondary organic aerosols and cloud condensation nuclei would increase regionally. Overall, we demonstrated that even a relatively mild scenario of tree species change will alter the energy balance and air quality in a way that could potentially influence the climate on a landscape scale.

1. Introduction

Climate change will have multiple impacts on temperate forests’ monospecific stands, which will be particularly vulnerable [1,2,3]. Therefore, forest management aims to increase the proportion of mixed forests in order to improve the resilience of future forests [4,5]. Based on this expectation, it is essential to investigate various stand structures and species compositions, and to evaluate their growth performance and susceptibility to future stressors [2,6,7,8]. Changes in forest structure and species composition, however, interact with the regional environment and affect the radiation budget, water balance through cloud condensation nuclei [9], and air quality through deposition of pollutants and emission of compounds [10]. While these feedbacks are known to be crucial on the regional to global scale [11,12,13], and the processes are often derived from studies at the local scale, the impacts that can be expected also on subregional scales due to altered species composition are not well investigated.

Introducing Silver fir (Abies alba Mill.) into European beech (Fagus sylvatica L.) stands has been proposed as a promising strategy to increase resilience of forests in Central Europe [14]. Silver fir, indicated as fir in the following, is less sensitive to spring drought events [15] and is generally expected to thrive under a warmer climate with somewhat less precipitation [16]. Thus, it can be expected that more wood growth and thus more carbon sequestration will be possible with than without implementing increased fir abundance, with beneficial impacts on general greenhouse gas balance. In addition, fir admixture may well increase the soil organic carbon content relative to pure beech forests, further increasing the long-term sequestration of the greenhouse gas CO₂ [17]. On the other hand, an increase in the fraction of fir relative to European beech, indicated as beech in the following, may lead to impacts on air chemistry and the subregional climate (see Figure 1 for illustration), for example, due to the different albedos of deciduous and evergreen species. Also, annual evapotranspiration has been estimated to be slightly higher, resulting from deep rooting as well as a high interception [18,19]. A major impact on air chemistry may also result from changes in the emission of biogenic volatile organic compounds (BVOCs) such as isoprene (ISO), mono- (MT), and sesquiterpenes (SQT) since beech and fir have different emission patterns regarding MT and SQT [20,21]. While it has been demonstrated to be absent in beech e.g., [22], it is still a debated question for fir [21,23]. In addition, emission patterns may respond differently to stress [24,25]. While terpenes tend to increase aerosol formation and thus the abundance of cloud condensation nuclei, ISO emission generally favors the formation of ozone and may reduce aerosol mass yield and cloud condensation nuclei from terpenes [26,27,28].

In the present work, we studied the influence of intermixing fir into pure beech forests on radiation balances and resulting changes in air temperatures, evaporation, and emission regime. The investigation was based on several years of research at a number of monitoring plots in the Black Forest, Southwestern Germany, where Silver fir constitutes an autochthonous tree species. Furthermore, we aimed to evaluate the secondary effects of ozone and aerosol formation and discuss the feedbacks on the radiation budget, air quality, and vegetation, e.g., due to potentially dangerous ozone levels. While decisions on tree species’ selections have previously been based on growth, productivity, and susceptibility to disturbances [14,29,30], this study supports additional considerations of regional climatic feedback.

2. Materials and Methods

2.1. Investigated Area

The investigations were carried out in the Black Forest, which is situated in the Federal State of Baden-Württemberg, Southwestern Germany (Figure 2, right). The Black Forest holds various tree species with Norway spruce (Picea abies L.) being the most abundant and beech and fir contributing considerable fractions (

Table 1. Tree species’ information for the Black Forest used for model estimates on the regional scale. Species and regionally averaged foliage biomass originate from Kändler and Cullmann [33] and the Forest Research Agency, Freiburg (pers. com.). Albedo and emission parameters (standard emission rate ER per dry weight (DW) for isoprene, monoterpenes, and sesquiterpenes, as well as temperature response β) were taken from Bonn et al. [24] (except bold numbers that were determined in this study).

Tree Species	Area Fraction [%]	Foliage Dry Weight [g/m2]	Albedo	ERISO [µgC gDW−1 h−1]/β ISO [K−1]	ERMT [µgC gDW−1 h−1]/β MT [K−1]	ERSQT [µgC gDW−1 h−1]/β SQT [K−1]
Picea abies	42.8	177 ± 31	0.09	0.05	0.886/0.11	0.35/0.11
Abies alba	18.5	174 ± 54	0.076	0.04	22.24/0.14	0.13/0.04
Fagus sylvatica	15.3	180 ± 14	0.13	0	49.1/0.31	0.1/0.11
Pseudotsuga menziesii	5.1	141 ± 45	0.09	0.008	0.064/0.26	0/-
Pinus sylvestris	5.0	77 ± 39	0.101	0.1	1.4/0.09	0.28/0.188
Quercus robur/petraea	3.1	179 ± 83	0.18	62.8	13.1/0.121	0.6/0.11
Betula spp.	1.9	105 ± 31	0.167	0	0.7/0.09	0.31/0.18
Alnus spp.	1.2	109 ± 92	0.16	0.018	0.13/0.09	0.0044/0.145
Carpinus betulus	1.2	132 ± 110	0.17	0.1	0.0093/0.09	0.011/0.145
Tilia spp.	1.2	110 ± 92	0.17	0.077	0.51/0.09	0.0067/0.145
Acer spp.	1.2	132 ± 49	0.17	0.1	1.8/0.09	0/ -
Populus spp.	1.2	85 ± 31	0.18	67.3	3.45/0.09	0/ -

). Experimental field studies (sap flow and exchange of trace gases’ measurements) were carried out at Freiamt (48.15° N, 7.91° E, ca. 400 m above sea level (a.s.l.)). Two different plots within this area were chosen to investigate the impact of fir admixture as possible future stand composition [14,19,30]: One is a monospecific beech forest (BB), the second is a forest in which beech and fir contribute approximately equal to total basal trunk area (BF). A more detailed description of both sites is provided in

Table 2. Properties of the pure beech stand (BB) and the mixed beech-fir stand BF sites as measured in 2017 (Da: Arithmetically averaged tree diameter at breast height [cm]. Dg: Tree basal area based on average diameter at breast height [cm]. Ha: Arithmetically averaged tree height [m]. Hg: Tree height of Dg [m]. N: Number of trees. BA: Tree basal area in 1.3 m height [m²]. VOL: Total stem volume [m³. Wtot: Total woody aboveground biomass (tons)). The data were obtained on plots of 300 m². Hg was calculated with logarithmic best-fit functions based on available height and diameter measurements; VOL and Wtot were estimated from Hg and Dg following Gasparini et al. [34].

	Pure Stand (BB)	Mixed Stand (BF)
	Beech	Beech	Fir	Total
Da	24.3			22.9
Dg	25.4	18.6	32.1
Ha	21.2			18.1
Hg	22.9	21.4	21.4
N ha−1	567	400	200	600
BA ha−1	28.6	10.8	19.9	30.7
VOL ha−1	261	77	170	247
Wtot ha−1	159	48	91	139

. Surface albedo of fir was determined at a rather pure fir stand at Wittnau (47.93° N, 7.84° E, ca. 650 m a.s.l.).

Mean annual air temperature and total annual amount of precipitation are 9.6 ± 6.5 °C and 1100 ± 554 mm (1984–2014, station 1440, Freiamt-Ottoschwanden, German Weather Service), respectively. The soil has a sandy loam structure that originates from a red sandstone and was classified as a Dystric Cambisol (according to the Food and Agriculture Organization of the United Nations (FAO) Soil Unit Classification [31]) with sandy loam texture [17]. A climate station for long-term observations set up by MeteoGroup (Berlin, Germany) in 2015 is situated at 440 m a.s.l., about 700 m north of the site at a cleared area. It supplied basic meteorological data in 20-min time resolution. Annual mean global radiation amounted to 135 ± 13 W m⁻², mean annual air temperature to 10.8 ± 0.8 °C, and total annual mean precipitation to 1010 ± 17 mm yr⁻¹ (2016: 997 mm yr⁻¹, 2017: 1022 mm yr⁻¹). For occasional gap filling and quality check, three nearby German Weather Service sites were used (Freiburg, Emmendingen-Mundingen, and Elzach-Fisnacht, https://www.dwd.de/DE/klimaumwelt/cdc/cdc_node.html). Additional information on air pollutants, particularly concentrations of O₃ and nitrogen oxides (NO_x), were obtained from the State Institute of the regional Baden-Württemberg office for Environment (LUBW) at the station Schwarzwald-Süd (47.81° N, 7.76° E, 905 m a.s.l.). Additionally, aerosol particle characteristics were obtained from the German Ultrafine Aerosol Network GUAN, [32]]) station at Schauinsland (47.91°N, 7.91°E, 1205 m a.s.l.) (for all station locations see Figure 2, left).

2.2. Radiation Budget Overview

The radiation budget at the top of the forest was calculated based on the following considerations: Subsequent to atmospheric interactions, a fraction of solar radiation reaches the Earth’s surface as global radiation (F_global), part of which is reflected (albedo α), while the rest (1 − α) is absorbed. The energy absorbed is converted into latent heat (F_ET, phase transfer of water), sensible heat (F_H), and into changes of the heat storage within the vegetation (dH_veg) and the soil (dH_soil). Only a small amount of F_global (<2% of net radiation budget) is directly used by photosynthesis [35], which was neglected henceforth. Therefore, the radiation budget on the local scale can be written as:

$$F_{global}·\left(1-\alpha \right)=F_H+F_{ET}+d{H}_{soil}+d{H}_{veg}$$

(1)

At the regional scale, the input of energy is better described by differentiating direct global radiation (F_{global_d}) and the energy that originates from longwave downward radiation (F_atmos) that has been absorbed and re-emitted by greenhouse gases, clouds, and aerosol particles in the atmosphere.

$$F_{global\_d}·\left(1-\alpha \right)+F_{atmos}=F_H+F_{ET}+d{H}_{soil}+d{H}_{veg}$$

(2)

With

$$F_{atmos}=\left(1+0.26·N\right)·(1-\left(1+\omega \right)·\exp (-{\left(1.2+3·\omega \right)}^{0.5}))·\sigma ·T_a^4$$

(3)

considering the humidity correction factor ω according to Prata [36]

$$\omega =46.5·\frac{vp}{T_a}$$

(4)

with vp being water vapor pressure (hectopascal), T_a ambient temperature (K), and σ the Stefan–Boltzmann constant (5.67·× 10⁻⁸ W m⁻² K⁻⁴). N is the cloudiness-covered fractions from zero to eight.

With F_{global_d} defined, the terms albedo α, latent heat F_ET, heat storage change (dH_veg, dH_soil), and sensible heat (F_H) need to be calculated in dependence on vegetation properties and, thus, on species’ composition. In order to solve the radiation budget, we newly measured α for firs and used literature data regarding beeches (Section 2.3). Although heat storage change can be regarded as negligible in the long term, it has been calculated to consider changes in seasonality based on directly measured differences in soil temperature and biomass storage (Section 2.4). F_ET is determined from mechanistic models that were evaluated with sap-flow measurements (Section 2.5). In addition, the effect of species’ composition change on F_atmos was estimated, considering the change in cloud properties and ozone concentration that can be expected with different BVOC emissions (Section 2.6). Therefore, the sensible heat flux F_H can be calculated iteratively, rearranging Equation (2):

$$F_H=F_{global\_d}·\left(1-\alpha \right)+F_{atmos}-F_{ET}-dH$$

(5)

As F_H is directly related to canopy surface temperature, a change in radiation flux can be transferred to a change in local temperature.

$$\Delta T_s={\left(\frac{\Delta F_{H,mix}}{\epsilon ·\sigma }\right)}^{\frac{1}{4}}-{\left(\frac{\Delta F_{H,beech}}{\epsilon ·\sigma }\right)}^{\frac{1}{4}}$$

(6)

with the emissivity ε assumed to be specific for deciduous trees (0.975) and conifers (0.99) [37]

2.3. Albedo

Surface albedo α of fir was derived from measurements of global radiation (F_global) and reflected shortwave radiation (F_reflect) using a component net radiometer (CNR4, Kippen & Zonen B.V., Delft, The Netherlands). The instrument was positioned at a forest plot in the vicinity of Wittnau (see Figure 2), which represented a closed canopy of pure firs that were about 30 years old and 14 m in height. The site had a negligible slope and was homogeneous across several hundred square meters, sufficient to avoid interferences by other species. We installed a mobile tower at the site with a platform approximately 1 m above the top of the stand. Radiation measurements were carried out between 28 July and 16 August 2017 with a time resolution of 60 s. Reflected radiation values smaller than 2 W m⁻² (global radiation of less than about 20–40 W m⁻²) were discarded from further analysis because of twilight conditions, which resulted in unrealistic α values. Albedo of beech was taken from literature (see Table 1).

2.4. Soil Temperature and Heat Storage Terms

Measurements of soil temperature and volumetric water content were carried out in 10, 15, 25, 50, and 80 cm depth at two-hour intervals starting in March 2016 until the end of November 2017 (applying Decagon 5TM VWC sensors). We used two sets of sensors that were implemented at the pure beech (BB) plot and two that were installed in the mixed plot in equal distance to beech and fir (BF). In addition, soil texture and other properties were analyzed from nine profiles throughout 1 m depth from which an average soil density of 1300 kg m⁻³ and a heat storage capacity of 2040 J kg⁻¹ K⁻¹ at water content at field capacity was determined (with no significant differences between plots). Water content was considered as affecting heat capacity and used along with the temperature differences to calculate heat storage changes throughout the year down to 1 m depth. For further information, also see Magh et al. [18]

In order to calculate aboveground heat storage, biomass of both stand types was estimated based on tree density and dimension measured at the test sites (see Table 2). Since the share of foliage biomass is highly variable during the season but usually less than 2%, aboveground biomass was assumed to be equal to woody biomass, and specific heat capacity was assumed as 2858 J K⁻¹ kg⁻¹ for beech and 3546 J K⁻¹ kg⁻¹ for fir (derived from the specific heat of water (4186 J K⁻¹ kg⁻¹), dry wood for both species (2021 and 2720 J K⁻¹ kg⁻¹, https://www.schweizer-fn.de/stoff/wkapazitaet/wkapazitaet_baustoff_erde.php, visited at 18 August 2020), and the assumption that average moisture content of woody biomass was 63 and 129%, for beech and fir, respectively [38]. Finally, the heat storage change was calculated considering local air temperature measurements after Thom [39].

2.5. Evapotranspiration and Latent Heat Flux

Latent heat (F_ET) was calculated from evapotranspiration by multiplication with the specific evaporation heat L at the trees’ surface that was determined separately for beech and fir. The results were averaged for the BF site. Overall, three model approaches were applied which were (1) the Penman–Monteith method (PMM) directly [40,41,42], (2) the model ‘Growth Of Trees Is Limited by Water’ (GOTILWA+) [43], which uses the Penman–Monteith formula for determining potential evapotranspiration but accounts for phenology as well as limitations in water supply, and (3) the Landscape-DeNitrification and DeComposition model (LDNDC) [44], which calculates potential evapotranspiration using a modified Thornthwaite approach [45,46,47]. As part of the water balance, GOTILWA+ and LDNDC also calculate transpiration, which are compared with measurements of sap flow carried out at the Freiamt site with the heat ratio method [48] for selected trees (three beech and two fir trees measured during July 2016). For more information about the methodology, see Magh et al. [18]. For better comparability to other terms, the evaluation of transpiration fluxes was used to check for the reliability of water budgets in order to use the PMM approach for the determination of F_ET in 2017.

2.6. Secondary Radiative Effects

2.6.1. Emission of Biogenic Volatile Organic Compounds

The emission of BVOCs from fir and beech were measured campaign-wise for 10 days in 2016 and 2017 during summer (June–July) and autumn (September–October) at the Freiamt site. Leaf and needle emissions were determined for sunlit branches at the upper canopy height (approx. 20 m) using a cherry picker and a GFS-3000 system (Walz, Germany), as described in Kreuzwieser et al. [49]. Cartridge samples were collected using air from sunlit branches enclosed in Nalophan foil for 60 min and of air samples nearby as references. BVOC concentrations were analyzed by GC-MS (see Yáñez-Serrano et al. [50] for details). Gas exchange was defined as the difference in concentration between enclosed branch air and reference samples obtained over enclosure time and per dry weight of green biomass. Results were used to define standard emission rates ER per dry weight (DW) and response coefficients β (K⁻¹) (according to Guenther et al. [51]) in order to calculate BVOC emission as a function of surface temperature T and photosynthetic active radiation I (ISO) or only T (MT, SQT).

$$E_{ISO}=E{R}_{ISO}\times C_{LT}$$

(7)

$$E_{MT}=E{R}_{MTS}\times exp\left(\beta \left(\mathrm{T}-T_S\right)\right)+E{R}_{MT}\times C_{LT}$$

(8)

$$C_{LT}=\frac{exp\left(c{T}_1\frac{\left(\mathrm{T}-T_S\right)}{R{T}_{S}T}\right)}{1+exp\left(c{T}_2\frac{\left(\mathrm{T}-T_M\right)}{R{T}_{S}T}\right)}\times \frac{\alpha c{L}_1\mathrm{I}}{{\left(1+{\alpha }^2{\mathrm{I}}^2\right)}^{0.5}}\times {\gamma }_{SM}$$

(9)

with empirical constants cT₁ = 95,000, cT₂ = 230,000, α = 0.0027, and cL₁ = 1.066, and T_M = 314 (K), as well as the standard temperature defined as T_S = 303 K. R is the universal gas constant (8.314 J K⁻¹ mol⁻¹) and T and I are leaf surface temperature in K and absorbed photosynthetic active radiation in µmol m⁻² s⁻¹, respectively. Other compounds than isoprene and monoterpenes were calculated assuming no light dependence at all. The factor γ_SM describes the dependency on soil moisture, which is a linear decrease with decreasing relative soil moisture below a certain threshold [25,52]. For the regional emission estimates, standard emission rates for other species in the area were taken from literature and scaled with area fraction and the dry mass density of leaves (Table 1).

2.6.2. Ozone and Secondary Organic Aerosol (SOA) Production

With BVOC emission parametrized using the measured data from the campaigns, the amount of O₃ production and the abundant O₃ volume mixing ratio were calculated with the box model ‘Chemistry As A Boxmodel Application’ combined with the ‘Module Efficiently Calculating the Chemistry of the Atmosphere’ CAABA/MECCA [53], which was extended by a BVOC emission and SOA scheme [54]. The applied BVOC emission scheme accounted for emissions of ISO, MT (assumed as a mixture of α-pinene, β-pinene, and limonene), and SQT (with the properties of β-caryophyllene). The basic equation for the formation of secondary organic compounds (J_ORG) is based on the description of 3-nm-sized particles in Bonn et al. [26]:

$$\begin{gathered} J_{ORG}=\left[SOZ\right]\times min\left(c{J}_1\left[RO2\right], c{J}_2\left[sCI\right]\right)\times c{J}_1\left[RO2\right] \\ +\frac{c{J}_2\left[sCI\right]}{kHO2\left[HO2\right]+kNO\left[NO\right]+kNO3\left[NO3\right]+kRO2\left[RO2\right]} \end{gathered}$$

(10)

$$SOZ=ts\times \frac{kOS\left[O3\right]\left[SQT\right]Y+kH2O\left[H2O\right]}{kOH\left[OH\right]+CS}$$

(11)

with SOZ describing the concentration of secondary ozonoids (particle-forming compounds) and RO2, sCI, HO2, NO, as well as NO3, are organic peroxy radicals, stabilized Criegee intermediates, hydroperoxyl radicals, nitrogen monoxide, and the nitrate radical, respectively, all derived from measurements and/or VOC emission simulations and used as concentrations in molecules cm⁻³. Reaction constants with units cm³ molecule⁻¹ s⁻¹ are all indicated as k_i,j and are taken from the Master Chemical Mechanism [55], as referenced in Bonn et al. [54]. Y denotes the stabilizing rate and CS the condensation sink of SOZ. The empirical constants are cJ₁ = 2.5 × 10⁻¹³ cm³ molec⁻¹ s⁻¹, cJ₂ = 2 × 10⁻¹² cm³ molec⁻¹ s⁻¹, and ts = 80 s⁻¹.

It was assumed that these emissions instantaneously mix within the meteorological boundary layer. Other components that were considered for chemical reactions are NO, NO₂, CO, benzene, and toluene concentrations (station Schwarzwald-Süd) as well as total particle mass and volatile fraction (assumed as predominantly organic) (station Schauinsland), provided in hourly time resolution. Background concentrations of further VOCs were obtained from additional proton-transfer-reaction mass spectrometry (PTR-MS) measurements (station Schauinsland, data from 25 July to 9 September 2017). Meteorological data of Freiamt station were used as input parameters (see Section 2.1). Simulations targeting tropospheric O₃ concentration, particle formation rates up to 3 nm in diameter, and SOA mass were carried out for the years 2016 and 2017, assuming (1) the current tree species’ mixture of the Black Forest according to Table 1, (2) medium increase of fir abundance, replacing one-third of the area currently occupied with beech, and (3) high increase of fir by replacing half of the current areas currently occupied by beech.

2.6.3. Cloud Condensation Nuclei and Cloud Properties

Subsequently, we estimated the impact of tropospheric ozone and SOA concentration on cloud properties. The effect of cloud properties was calculated according to Spracklen et al. [56]:

$$\mathsf{\Delta }{F}_c=-F_{solar}·A_c·T{R}^2·\mathsf{\Delta }{R}_c$$

(12)

ΔF_c describes the change of all radiation fluxes from clouds that is induced by the change in cloud condensation nuclei CCN. A_c represents the fractional cloud cover, the median of which is ${0.75}_{0.13}^{0.88}$ (Mt. Feldberg, German Weather Service, station no. 1443). TR abbreviates the fractional atmospheric transmission of shortwave radiation above the cloud layer, assumed as 0.76 [56], and ΔR_c is the change in cloud albedo, calculated as

$$\Delta R_c =\frac{1}{3}{R}_c·\left(1-R_c\right)·\frac{\mathsf{\Delta }\mathrm{CCN}}{CCN}$$

(13)

The cloud albedo R_c was assumed as 0.20 ± 0.05 [57]. CCN is determined using an empirical relationship with the emission of monoterpenes (ER_MT) that was suggested previously based on measurements in a very similar structured area [58]:

CCN (s = 0.3%) = 675 ± 31 cm⁻³ + 0.45 ± 0.10 g(DW) h µg⁻¹ cm⁻³ ER_MT

(14)

(This equation is limited in applicability for CCN concentrations less than ca. 800 cm⁻³, which is, however, well above 97% of all measurements at the investigated site (station Schauinsland)).

3. Results and Discussion

3.1. Land Surface Albedo

The measurement period covered a wide range of different meteorological conditions with 31,600 data points in total. Plotting surface albedo vs. incoming global radiation provided an equilibrium value of 0.076 ± 0.001 at sufficiently high radiation (Figure 3). This value is substantially lower than the literature-based albedo value of beech of 0.13 ± 0.01 [59,60]. In addition, leaf area of beech is negligible between early November and late April (with a peak leaf area index of 6 and transition phases of about four weeks in spring and autumn), while the evergreen fir provides continuous foliage cover (leaf area index of 3.8 in winter and up to 9.8 in summer). These phenological results are comparable with data available for Central European conditions [61,62,63]. Overall, a pure canopy of fir thus yields a surplus of 7.3 W m⁻² as an annual mean and 54.4 W m⁻² at maximum compared with a pure beech canopy. Mixtures are to be found in between. In other words, fir increases the heat flow into the ecosystem, although not necessarily the temperature of the canopy (see next section).

3.2. Heat Storage Dynamics

Although soil temperature dynamics are usually less expressed than air temperature changes, the comparison between the BF and BB site shows that fir contribution results in a faster heating and cooling during the day. This is probably due to a faster energy transport enabled by the crown shape of firs (see inserted figure within Figure 4, top left). On a seasonal scale, soil at fir-admixed plots is cooler in winter, which leads to slightly less heat storage compared to pure beech forests of about 0.3 K in annual average (Figure 4, bottom left). The effect is accelerated by the finding that the soil under beech trees is generally wetter (see next section), leading to more heat storage in the ground of the BB site than under the BF site throughout most of the year (Figure 4, bottom right). The aboveground biomass, however, is counteracting this development, since the heat capacity of fir wood is larger. Although the overall biomass at the BF site is slightly smaller, which is due to the relatively many but small beech trees, this effect is still sufficient to result in up to 4 W m⁻² smaller heat storage in mid-summer at the BB site (Figure 4, bottom right). These findings show (1) that the increased heat flow into forests with fir admixture is quickly transported throughout the canopy to heat up the ground and not the foliage surface and (2) that not the ground but the trunks are the major sink of energy that results from a smaller albedo.

3.3. Latent Heat Exchange

Simulated transpiration for beech is characterized by the restriction to the growing season in both GOTILWA+ and LDNDC and by a general dependency on temperature (Figure 5). Evaluation of simulations by sap-flow measurements indicates that both models seem to be able to capture the response for beech, while GOTILWA+ overestimates transpiration for fir (Figure 6, left). Evapotranspiration, as represented by the PMM approach, is generally a bit larger for beech but smaller for fir than the fluxes derived by the two ecosystem models (Figure 6, right).

Similar results can be derived from the annual water flux budget. Evapotranspiration for the years 2016/2017 was calculated to be 574/627 (PMM), 407/504 (LDNDC), and 483/444 mm (GOTILWA+) for beech stands and 624/686 (PMM), 465/572 (LDNDC), and 479/441 mm (GOTILWA+) for fir. This corresponds to 59/65 (PMM), 45/51 (LDNDC), and 46/46 % (GOTILWA+) of the annual precipitation for beech and fir, respectively. For a watershed within the Black Forest, overall evapotranspiration rates have been determined to be 55% of precipitation, estimating the value to be 63% for densely covered forest patches [64]. Compared to these investigations, GOTILWA+ and LDNDC seem to slightly underestimate evapotranspiration while the PMM provides values in a reasonable range. Therefore, we used PMM, which can also easily be used for smaller time steps, for energy, and air chemistry calculations. However, we considered the results from both ecosystem models for determining the range of uncertainty.

Assuming that mixtures of beech and fir can be represented by patches of pure stands that are proportional to the basal-area contribution of the respective species within the mixed forest, evapotranspiration fluxes can be estimated for mixed stands. Thus, the 50% mixture of the BF site increases evapotranspiration by 59 mm per year, although it can even be reduced during the summertime. This corresponds to an annual median F_ET of −66 W m⁻². In other words, annual ET would rise and F_ET would decline by about 7.6% per 10% increase in fir intermixture. This has numerous implications. On the one hand, it partly counteracts the increased heat inflow due to the darker canopy, decreasing the effect on regional energy balance. On the other hand, however, increased water flow back into the air directly feeds back on incoming radiation directly and indirectly by affecting cloud formation.

3.4. Secondary Radiation Effects

3.4.1. Biogenic Volatile Organic Compounds

Measurements indicated that neither beech nor fir emitted significant amounts of isoprene, thereby confirming studies on beech [20,22] and fir in the French Vosges forest [21], which are in close vicinity to the Black Forest. However, our data contradict the rather high emissions reported for fir from Belgium [23]. For the whole area, ISO inputs into the box model calculations thus mainly originate from small contributions of high emitting oak and poplar trees at lower altitudes (Table 1). In contrast, emission patterns of MT and SQT were found to be considerably different between beech and fir, resulting in substantial variations in total emissions and reactivity: Emission rates for SQT under standard conditions (ER_SQT) were 0.2 µg g(DW) h⁻¹ for beech and 0.7 µg g(DW) h⁻¹ for fir. Temperature response coefficient β [51] of fir were, however, smaller (β_SQT^fir = 0.04 K⁻¹) than those for beech (β_SQT^beech = 0.11 K⁻¹), resulting in particularly large differences at low temperatures and convergence of SQT emissions when T_S approaches ≈ 45 °C. Together with the higher leaf area of fir, the overall response was an increasing emission of MT and SQTs with increasing fraction of firs throughout the year. For our sites, the actual SQT emission was calculated to be about 75% and the MT emissions 55% higher for intermixtures with 50% fir than for pure beech sites. On a regional scale, these differences led to a considerable increase of BVOC emissions with impacts on ozone and aerosol formation, as discussed in the next sections.

3.4.2. Ozone Formation

Tropospheric O₃ concentration [O₃] was calculated based on BVOC emissions and background measurements of NO_x and meteorology [65]. It was found that [O₃] in the western part of the Black Forest depends to a similar degree on NO (NO_x) and HO₂ (BVOC). As tropospheric O₃ is changed on a regional scale, we used either the current beech/fir mixture (current) or a scenario with a higher share of fir (admixing scenario) at all sites where beech forests are currently observed in the Black Forest, while the contribution of species other than beech or fir remained the same. Simulations based on current conditions resulted in O₃ mixing ratios that agreed fairly well with measurements (r² = 0.83) [54]. The overall radiation effect of mean tropospheric O₃ (F_O3) calculated for 2017 varied between 0 and 2 W m⁻² (Figure 7). The increase of tropospheric O₃ with increasing contribution of fir (BF scenario) resulted in a radiative change of up to +0.25 W m⁻² on an hourly basis. The mean annual increase of F_O3 was estimated to be moderate with 0.023 W m⁻².

An increase of fir abundance according to the scenario was thus simulated to cause an increase of mean tropospheric O₃ concentrations by +11.1%. Since the difference in emission depends on temperature (see previous section), the increase in [O₃] is smaller under cool conditions and larger when the weather is hot.

3.4.3. Secondary Organic Aerosols and Cloud Properties

The annual median of the direct radiation effect of SOA (F_SOA) for current conditions was calculated to be −0.089 W m⁻². This value was slightly less negative (−0.087 W m⁻²) for the tested scenario with increased fir abundance (Figure 8). While this effect is small and the degree of mixture is negligible, at highly polluted episodes with very warm days it can increase to almost −4 W m⁻². During these times the increased mixture scenario differs by more than 1 W m⁻², modifying the SOA-related radiation effect by 50%.

Further changes were approximated for formation rates of particles (up to 3 nm in diameter) with subsequent changes in concentrations of cloud condensation nuclei (CCN). This effect modifies the radiation balance of clouds, i.e., the brighter the clouds, the less radiation reaches the Earth’s surface. The reduction is a function of water vapor supersaturation necessary for CCN activation and also depends on monoterpene concentrations [58]. The radiation effects of CCN activation on the energy balance when modifying the forest composition are displayed in Figure 9 (also demonstrating the effect of water supersaturation). It shows that the effect was much larger than could be derived from SOA formation and that annual means are not very indicative. Since CCN formation underlies a pronounced seasonal cycle, the impact of fir admixture on cloud formation is especially expressed in summer.

3.5. Total Radiation Budget

The changes in radiation fluxes described in Section 3.1, Section 3.2, Section 3.3 and Section 3.4 were combined to estimate the average annual impact of a change in species’ composition on the energy budget at the upper forest canopy. An overview of mean values of individual effects and strengths is listed in

Table 3. Mean annual radiation contributions based on the year 2017 at Black Forest for the current tree distribution and an increased fir contribution.

Radiation Effect	Scale *	BB Site (L) or Present Mixture (R) [W m−2]	Intermixture of 50% Fir into All Current Beech Stands (BF, Scenario) [W m−2]	Difference between Present Mixture and Scenario [W m−2]
Fglobal_d·(1 − α)	L	105.6	112.2	6.6
FET	L	−48.0	−66.1	−18.1
Fatmos	R	299.8	299.8	0.000
FO3 **	R	0.284	0.308	0.023
FSOA **	R	−0.089	−0.087	0.002
FCCN **	R	0.769	0.749	−0.021
ΣFi = −FH ***	L	357.4	345.8	−11.6
Tmean (°C) ****	L	10.4	7.0	−3.4

* L—local, R—regional. ** Radiation fluxes, which are part of the energy that originates from longwave downward radiation (F_atmos) for total budgeting. *** Sum of all energy fluxes. **** Canopy temperature.

. The values obtained are similar in magnitude to Holst and Mayer [66] for a nearby beech site and to Schmid et al. [67] for a North American deciduous forest. While net shortwave radiation input at the canopy top tentatively increased with increased intermixture of fir trees from 106 to 112 W m⁻², F_ET tended to decrease from −48 to −66 W m⁻². However, a high uncertainty of more than 100% in F_ET resulted from the large variability of estimated evapotranspiration. The overall heat stored in the ecosystem increased with increasing fir admixture, although the uppermost canopy was simulated to even cool down. While the differences in radiation budget between the current conditions and the scenario were relatively small on an annual scale, species’ mixture effects can be more expressed at shorter timescales, i.e., days. This means that our research supports the assumption that species’ mixture may become particularly important during extreme heat periods, leading to intensified short-term stresses that are often more damaging than cumulative ones [68].

Large uncertainties derive from daily climate variations, for instance incoming global radiation that triggers specific components of the radiation balance differently [69,70]. As demonstrated here for summer periods, also Oliphant et al. [70] measured substantial radiative contributions of dH_veg for a deciduous forest at daily and sub-daily scale with up to 10 W m⁻², resulting in crown surface temperate increases of up to 4 K at noon. Similarly, upper canopy surface temperatures of beech trees were found to be 4–5 K higher than air temperature, while canopy and air temperatures of conifers were approximately the same [71]. Other results confirm that conifer forests considerably increase their heat storage during summer [35]. It should be noted that these differences may represent a feedback response to BVOC emission, modifying the emission conditions differently in different parts of the canopy, which was, however, not considered here.

3.6. Uncertainties and Deficiencies of the Investigations

First, we want to emphasize that the investigation did not consider any disturbances such as damages by herbivory and browsing, which might penalize Silver fir compared to competing tree species. For example, it has been shown that firs were preferentially browsed in comparison to Norway spruce [72], although this also depends on the forest mixture and the abundance of other species within the understory [72,73]. Second, processes that were considered as, for example, describing the water balance, may include aspects that could not be considered here, such as deep water redistribution [74]. This may compensate for the higher evaporation losses from interception that results from the higher potential leaf area of firs. Even if considered, as in the case of the drought impact on BVOC emission, part of the response such as an increase of emission under mild drought is still lacking a consistent mathematical description [75]. Third, many of our investigations were based on local measurements that we assumed to be representative for the applied region in order to present a consistent scaling investigation. However, only few investigation sites of limited size were employed, which, together with not differentiating between biomass in wood and foliage, resulted in a relatively high uncertainty, in particular, regarding the biomass at the sites. It would, thus, be favorable if more fir sites could be incorporated into a future scaling exercise. Finally, it should be noted that model results were still far from perfectly evaluated. We tried to compensate for this problem by using several approaches that, however, still might need improvements, e.g., regarding soil and forest structural considerations [76,77].

4. Conclusions

According to our results, a scenario assuming a 50% admixture of fir into all existing beech forests would lead to higher shortwave radiation input because of the darker foliage. The surplus energy is stored dominantly in the trunks of the conifers while the uppermost canopy layers may even cool due to the increased turbulence of the canopy space of firs. On a diurnal basis, this effect also leads to a relative cooling at daytime but warming during the night. The other important effect is the loss of energy and water by increased evapotranspiration, predominantly driven by a higher interception. This interacts with energy terms since less water decreases the heat storage capacity of the soil.

On the regional scale, MT and SQT emissions were simulated to increase by around 10%, resulting from 55–75% higher emissions at the local scale. This was simulated to increase O₃ as well as SOA production with implications for CCN concentration. Together with increased water evaporation, CCN leads to qualitative and quantitative changes of clouds, i.e., cloud brightness. Both air chemistry impacts and cloud formation effects have, in turn, impacts on the regional radiation regime and rainfall pattern and, thus, represent a climate feedback process, which, however, was not further evaluated in this study.

Despite the relatively minor change that was assumed regarding the scenario species’ composition, we overall estimated primary as well as secondary impacts of a species’ change on regional radiation balance and air chemistry to be substantial. In this particular case, increased energy absorption and evapotranspiration may increase the drought stress and need to be weighed against a possible higher resistance that fir may have compared to beech [14,15]. On the regional scale, an increase in [O₃] may be generally negative but tolerable or not depending on the specific conditions, while an increase in cloudiness, brightness of clouds, or even precipitation could be regarded as a positive climatic feedback or as an effect that decreases radiation availability for vegetation productivity. This, however, can only be judged in long-term, coupled-model investigations that consider these feedbacks.