Marc Peaucelle
email: marc.peaucelle@inrae.fr

Josep Peñuelas
email: josep.penuelas@uab.cat

Hans Verbeeck
email: hans.verbeeck@ugent.be

Accurate phenology analyses

Keywords: plant phenology, budburst, temperature, light, energy budget, modelling, climate warming

HAL is

Introduction

Plant phenology, the study of the timing of life-cycle events, drives several ecosystem functions, such as plant productivity and biomass, but also local and global climates by affecting biogeochemical and biogeophysical processes, such as carbon storage and energy fluxes [START_REF] Peñuelas | Responses to a warming world[END_REF][START_REF] Peñuelas | Phenology feedbacks on climate change[END_REF] , and the abundance and diversity of local flora and fauna, such as pollinators and herbivores [START_REF] Ramos-Jiliberto | Phenology determines the robustness of plant-pollinator networks[END_REF][START_REF] Chuine | Why does phenology drive species distribution? Philosophical transactions of the Royal Society of London[END_REF] . Understanding the environmental controls and responses of plant phenology to climate change is thus essential for several sectors, e.g. agriculture, forestry and gardening [START_REF] Chmielewski | Phenology in Agriculture and Horticulture[END_REF] , but also for conservation [START_REF] Morellato | Linking plant phenology to conservation biology[END_REF] and public health [START_REF] Katelaris | Climate change: allergens and allergic diseases[END_REF] (e.g. allergies).

It is now largely assumed that bud-break is induced by warming air temperature during spring [START_REF]Phenology: An Integrative Environmental Science[END_REF] in temperate and boreal regions, and this is the reason why phenology assessments mainly rely on critical air temperature sums preceding leaf unfolding, often referred to as the Growing Degree Day concept.

Climatic warming has strongly shifted phenophases in the Northern Hemisphere in recent decades [START_REF] Peñuelas | Responses to a warming world[END_REF][START_REF] Cleland | Shifting plant phenology in response to global change[END_REF][START_REF] Fu | Recent spring phenology shifts in western Central Europe based on multiscale observations[END_REF][START_REF] Jeong | Phenology shifts at start vs. end of growing season in temperate vegetation over the Northern Hemisphere for the period 1982-2008[END_REF] .

Rising temperatures have lengthened the annual growth cycle by advancing leaf unfolding in spring and delaying leaf fall in autumn [START_REF] Liu | Delayed autumn phenology in the Northern Hemisphere is related to change in both climate and spring phenology[END_REF] , albeit with variations among species [START_REF] Vitasse | Leaf phenology sensitivity to temperature in European trees: Do within-species populations exhibit similar responses?[END_REF] and regions [START_REF] Wang | Temporal Trends and Spatial Variability of Vegetation Phenology over the Northern Hemisphere during 1982-2012[END_REF] . Recent evidence, though, suggests that the sensitivity of spring phenology to warming is decreasing in northern forests [START_REF] Fu | Declining global warming effects on the phenology of spring leaf unfolding[END_REF] and that the rate of change in plant productivity does not match that of air temperature [START_REF] Huang | Velocity of change in vegetation productivity over northern high latitudes[END_REF] . Indeed, plant phenology may be acclimated to long-term biogeographical constraints [START_REF] Peaucelle | Spatial variance of spring phenology in temperate deciduous forests is constrained by background climatic conditions[END_REF][START_REF] Zohner | Interactive climate factors restrict future increases in spring productivity of temperate and boreal trees[END_REF][START_REF] Montgomery | Phenological responses of temperate and boreal trees to warming depend on ambient spring temperatures, leaf habit, and geographic range[END_REF] and may be co-limited by several other factors, such as light [START_REF] Zohner | Day length unlikely to constrain climate-driven shifts in leaf-out times of northern woody plants[END_REF] , water [START_REF] Peñuelas | Complex spatiotemporal phenological shifts as a response to rainfall changes[END_REF][START_REF] Papagiannopoulou | Vegetation anomalies caused by antecedent precipitation in most of the world[END_REF][START_REF] Delpierre | Modelling interannual and spatial variability of leaf senescence for three deciduous tree species in France[END_REF] and nutrients [START_REF] Fu | Nutrient availability alters the correlation between spring leaf-out and autumn leaf senescence dates[END_REF] . These observations suggests that warming does not have the same effect everywhere [START_REF] Seyednasrollah | Leaf phenology paradox: Why warming matters most where it is already warm[END_REF] , which has increased interest in other environmental drivers in recent decades, especially illustrated by multiple debates about the specific role of light (and photoperiodism) in spring phenology [START_REF] Zohner | Day length unlikely to constrain climate-driven shifts in leaf-out times of northern woody plants[END_REF][START_REF] Chuine | Warming, photoperiods, and tree phenology[END_REF][START_REF] Vitasse | What role for photoperiod in the bud burst phenology of European beech[END_REF][START_REF] Way | Photoperiod constraints on tree phenology, performance and migration in a warming world[END_REF][START_REF] Caffarra | Modelling the timing of Betula pubescens budburst. II. Integrating complex effects of photoperiod into process-based models[END_REF][START_REF] Körner | Plant science. Phenology under global warming[END_REF][START_REF] Fu | Daylength helps temperate deciduous trees to leaf-out at the optimal time[END_REF][START_REF] Singh | Photoperiod-and temperature-mediated control of phenology in trees -a molecular perspective[END_REF][START_REF] Flynn | Temperature and photoperiod drive spring phenology across all species in a temperate forest community[END_REF] .

How light affects spring phenology remains an open question. Most commonly, its effect is considered via photoperiod, often referred to as daylength. The daylength hypothesis implies that the quality and/or quantity of light is somehow directly sensed by plants through biochemical mechanisms.

Some recent studies suggest that the spectral composition of light can indeed influence foliar phenology [START_REF] Brelsford | The influence of spectral composition on spring and autumn phenology in trees[END_REF][START_REF] Strømme | UV-B and temperature enhancement affect spring and autumn phenology in Populus tremula[END_REF] . Light also plays a key role in regulating phytohormones, but the underlying mechanisms remain unknown [START_REF] Singh | Photoperiod-and temperature-mediated control of phenology in trees -a molecular perspective[END_REF] and clearly require more investigation. More sporadically, the effect of light has been treated as the sum of insolation over a specific period [START_REF] Fu | Increased heat requirement for leaf flushing in temperate woody species over 1980-2012: effects of chilling, precipitation and insolation[END_REF] , for which plants need a specific quotum for a phenological event to occur. The quantity and quality of light depend on plant location, which is the main reason why a response to daylength has often been proposed as a safety mechanism against frost at high latitudes and elevations. Only 35% of the woody species in the Northern Hemisphere, however, depend on daylength as a direct signal for leaf-out [START_REF] Zohner | Day length unlikely to constrain climate-driven shifts in leaf-out times of northern woody plants[END_REF] , and these species are mainly at mid-to low latitudes.

Light effect on spring phenology is still being debated. Recent studies nonetheless suggest a complex interaction between temperature and light. Daytime and nighttime temperatures during winter and spring have an asymmetrical effect on leaf unfolding [START_REF] Huang | Effect of preseason diurnal temperature range on the start of vegetation growing season in the Northern Hemisphere[END_REF][START_REF] Meng | Opposite effects of winter day and night temperature changes on early phenophases[END_REF][START_REF] Zhang | Responses of bud-break phenology to daily-asymmetric warming: daytime warming intensifies the advancement of bud break[END_REF][START_REF] Meng | Divergent responses of spring phenology to daytime and nighttime warming[END_REF][START_REF] Bigler | Daily Maximum Temperatures Induce Lagged Effects on Leaf Unfolding in Temperate Woody Species Across Large Elevational Gradients[END_REF] , with a greater weight of temperature during the day [START_REF] Meng | Opposite effects of winter day and night temperature changes on early phenophases[END_REF][START_REF] Fu | Three times greater weight of daytime than of night-time temperature on leaf unfolding phenology in temperate trees[END_REF][START_REF] Piao | Leaf onset in the northern hemisphere triggered by daytime temperature[END_REF] . Whether or not plants are able to sense light, radiation has a physical impact on plants: it affects the temperatures of their tissues. Since temperature has been shown to be sensed at the bud level [START_REF] Vitasse | Impact of microclimatic conditions and resource availability on spring and autumn phenology of temperate tree seedlings[END_REF][START_REF] Azeez | EARLY BUD-BREAK 1 and EARLY BUD-BREAK 3 control resumption of poplar growth after winter dormancy[END_REF] , omitting the physical effect of radiation introduces large biases into the interpretation of spring phenological responses based on air temperature.

In their recent study, Vitasse et al. highlighted the strong phenological effect of bud albedo and light exposure, explaining shifts in the budburst date reaching up to 12 days [START_REF] Vitasse | Impact of microclimatic conditions and resource availability on spring and autumn phenology of temperate tree seedlings[END_REF] , and revealing an important role of microclimatic variation on phenology. Indeed, bud temperature (Tbud) depends on its energy balance [START_REF] Hamer | The heat balance of apple buds and blossoms. Part I. Heat transfer in the outdoor environment[END_REF] . During the day, plant tissues absorb both shortwave (SW, visible and near-infrared) and longwave (LW, infrared) radiation from the sky but also radiation emitted and reflected by the surrounding environment (vegetation, soil) (Fig. 1a). Only a fraction (α, absorptivity) of SW radiation will be absorbed depending on bud traits such as color, coating, shape and size (Fig. 1b), while most LW radiation will be absorbed by buds. According to the Stefan-Boltzmann law, buds lose energy via LW radiation emission, while they absorb LW radiation emitted from surrounding objects. Finally, a part of their energy is lost by conduction and mostly by convection [START_REF] Landsberg | Apple bud and blossom temperatures[END_REF][START_REF] Grace | The temperature of buds may be higher than you thought[END_REF] (e.g. due to wind) while leaves lose an important part of their energy via transpiration. Tbud increases when energy gains exceed losses (Fig. 2a) and vice-versa. Tbud can thus be lower than air temperatures (Tair) on clear nights [START_REF] Landsberg | Apple bud and blossom temperatures[END_REF] or because of wind.

On the other hand, Tbud can be significantly higher than Tair during the day. The link between Tbud and energy balance has been known for more than 30 years [START_REF] Hamer | The heat balance of apple buds and blossoms. Part I. Heat transfer in the outdoor environment[END_REF][START_REF] Landsberg | Apple bud and blossom temperatures[END_REF] . Since then, all major studies linking temperature and photoperiod to phenological changes, however, have not accounted for the true temperature of plant organs.

What can we expect if we account for micrometeorology and bud temperature in phenological studies? Unfortunately, the lack of in situ observations for bud temperature does not allow to answer directly this question. As part of the reflection, we thus applied existing well-established energy balance approaches [START_REF] Hamer | The heat balance of apple buds and blossoms. Part I. Heat transfer in the outdoor environment[END_REF][START_REF] Landsberg | Apple bud and blossom temperatures[END_REF][START_REF] Muir | tealeaves: an R package for modelling leaf temperature using energy budgets[END_REF] to explore the potential variability in temperature of an isolated bud. This situation is well representative of the conditions encountered by sun-exposed buds of a tree, and especially of deciduous species (i.e. with no or minimum shading). As commonly applied in 'big-leaf" models where an entire canopy is represented by a single 'big' leaf, discussing the microclimate effect on an isolated bud will help us to explore the variability in phenology we can expect at different spatiotemporal scales and between species.

Results

Non-linear response of plant tissue temperature to microclimate

As a first example, we looked at the variability in Tbud estimated from its energy balance and site meteorological observations for an European Beech forest [START_REF] Granier | Ten years of fluxes and stand growth in a young beech forest at Hesse, North-eastern France[END_REF] . On average, Tbud is expected to be higher than Tair during the preseason (~1°C in our example; Fig. 2b). Day and night Tbud are higher or lower than Tair by several degrees. The temperature of buds thus strongly depends on the diurnal radiative cycle and the spectral composition of the light (SW/LW radiation), echoing the observed asymmetrical effect of diurnal temperatures on leaf unfolding [START_REF] Meng | Opposite effects of winter day and night temperature changes on early phenophases[END_REF][START_REF] Fu | Three times greater weight of daytime than of night-time temperature on leaf unfolding phenology in temperate trees[END_REF][START_REF] Piao | Leaf onset in the northern hemisphere triggered by daytime temperature[END_REF] . Applied on four other sites, this approach leads to similar results despite differences in Tbud profiles induced by differences in radiation along a latitudinal gradient (Fig. 3). Spring phenology does not only respond to average preseason temperature, but mainly to the accumulated effect of temperature and its dynamics. It is often assumed that chilling and forcing temperature required for budburst are only effective over specific windows, generally between 0 and 5 °C and over 5°C, respectively. Daily bud temperature variability might thus be the most important factor influencing leaf unfolding, not necessarily its average temperature. We could expect that the difference in extremum temperature sensed by buds over the preceding months (ΔTmin and ΔTmax, Fig. 2b) will inevitably affect the apparent forcing and chilling requirement for leaf unfolding.

Accounting for the energy budget of buds for six common species across Europe (Extended Fig. 1) we also expect a stronger interannual variability in Tbud than Tair, as well as different temporal evolutions over the last decades (Fig. 4). In our example, buds are expected to warm faster or slower than air depending on location and species, with 13 % and 7 % of the sites exhibiting an increase and a decrease in ΔT over 1990-2015, respectively. These trends represent idealized sun-exposed conditions without site or species-specific calibration. Still, we observe that the heterogeneity in ΔT evolution results from a complex and non-linear response to the amount of absorbed radiation and convection processes (Extended Fig. 2). Because leaf unfolding is earlier in 2015 than in 1990, the average amount of absorbed radiation during the preseason slightly decreased over this period, while most of the interannual variability in ΔT is driven by conduction and convection (i.e. wind). The difference in air-bud temperature and their non-linear and non-proportional relationship suggests that our current interpretation of the apparent bud sensitivity to warming needs to account for the temperature sensed by the plant.

Response to warming relies on organ traits and microclimate

We illustrated the role of bud energy balance through an idealized and constant representation of buds and their environment for all sites and species. Larger spatial and temporal variations are expected due to the effect of topography, ground albedo (e.g. snow, understory), differences in bud traits (Fig. 1b) and micrometeorological conditions 53 that will affect plant tissues energy balance. By affecting the amount of radiation reaching the buds (Fig. 1a), varying ground albedo from 0.1 (~wet bare soil) to 0.9 (~snow) leads to a doubling in preseason ΔT (Fig. 5). Since ground albedo strongly varies in space and over the preseason (e.g. snow), we can expect substantial differences in the phenological signal at the regional scale induced by radiation, as already observed from leaf unfolding observations [START_REF] Peaucelle | Spatial variance of spring phenology in temperate deciduous forests is constrained by background climatic conditions[END_REF] .

Different bud colors or coating will also affect solar absorption of specific wavelengths, while shape and size will modify convection processes and the amount of intercepted radiation (Fig. 6), and hence, bud temperature. For example, Common Ash (Fraxinus excelsior) has black buds while sycamore (Acer pseudoplatanus) has green buds and mountain ash (Sorbus aucuparia) have dense white trichomes (i.e. hairs) on their surfaces. In our example, a difference in solar absorptivity of 0.3 leads to a doubling in ΔT (Fig. 2a, Fig. 3a). The differences in bud traits can thus partly account for the observed interspecific differences in heat requirement and apparent sensitivity to temperature. This suggests that the phenological response of plants to warming might be more species-specific than we thought, which should be accounted for in large scale studies.

Discussion

Despite its central role at the organ level [START_REF] Zellweger | Forest microclimate dynamics drive plant responses to warming[END_REF][START_REF] Bailey | A new three-dimensional energy balance model for complex plant canopy geometries: Model development and improved validation strategies[END_REF] , micrometeorology is rarely accounted for in phenology studies because rarely measured, or simply because it is impossible to account for its effect such as in remote sensing analysis or terrestrial biosphere modelling. Instead, phenology studies, either local or regional, often use meteorological and climate datasets with hourly to daily time resolutions. The use of a steady state energy balance is easily justified under such conditions since thermal time constants of tree buds varies between a few seconds to about ten minutes [START_REF] Michaletz | A heat transfer model of crown scorch in forest fires[END_REF] . Accounting for average preseason radiation and wind conditions might better explain the observed variability in plant phenology than air temperature alone. Here, we only explored spring Tbud variability in the case of sun-exposed buds with no shading. Accounting for the potential protecting effect of leaves or needles in evergreen species might substantially attenuate the effects of radiation and wind on intra-and bottom-canopy buds. The concomitant use of high-resolution microclimate data and transient energy budget models will be needed to quantify such effects.

Drivers of phenological events and light are virtually impossible to separate, because daylength and radiation are strongly correlated with the time of year. Accounting for organ energy balances is thus promising for separating the environmental drivers of phenology using a single approach and potentially for reconciliating the differences observed in the field. Applying existing modelling approaches in the context of sun-exposed buds suggested that air temperature might be an imprecise and biased predictor of bud temperature and more importantly of its variability over the months preceding leaf unfolding, which might introduce biases in the analysis of chilling and forcing requirement for budburst. However, we also showed that bud temperature results from a complex combination of several biotic and abiotic factors, and under certain conditions air temperature might remain a good proxy for bud temperature.

The examples we have presented open new avenues to investigate and refine our interpretation of plant phenological acclimation to warming that were solely based on air temperature. For example, changes in spring radiation regimes over the last decades [START_REF] Sanchez-Lorenzo | Reassessment and update of long-term trends in downward surface shortwave radiation over Europe (1939-2012)[END_REF][START_REF] Pfeifroth | Trends and Variability of Surface Solar Radiation in Europe Based On Surface-and Satellite-Based Data Records[END_REF] have been hypothesized to increase discrepancies between standard air temperature and bud temperature [START_REF] Vitasse | Impact of microclimatic conditions and resource availability on spring and autumn phenology of temperate tree seedlings[END_REF] , which is in line with the energy budget theory described in this perspective. However, current observations do not allow such reassessment. Bud traits and in situ temperature observations are scarcely described in the literature. New experiments and observations are clearly needed for accurately quantifying the traits and energy budget of buds. Existing studies have mostly focused on leaves, but other organs should also be investigated. Key traits that will need to be measured to assess the interspecific variability of phenology include organ traits influencing solar absorptivity and heat storage, but organ temperatures (i.e. using thermocouples) concomitant with micrometeorological variables will also need to be directly measured. With the assumption that bud transpiration is negligible, their energy budget is simpler than for leaves. Properly calibrated, accounting for bud energy balance could improve the accuracy of phenological models that are still unable to predict the spatiotemporal variability of plant dynamics with satisfactory accuracy [START_REF] Richardson | Terrestrial biosphere models need better representation of vegetation phenology: results from the North American Carbon Program Site Synthesis[END_REF] . Such approach could not only improve spatiotemporal assessments of phenology (e.g. based on ground-based and remote sensing observations), but also provide new insights into ecosystem functioning such as the influence of global warming synchrony between flowers and insects, or between plant functions and climate (e.g. growth and soil moisture, signaling, etc.).

Finally, we stress that energy balance affects the temperature extrema sensed by plants (Fig. 2b).

The lengthening of the growing season in recent decades has also been associated with an increase in environmental risks. For example, earlier leaf unfolding exposes plants to late frost [START_REF] Liu | Extension of the growing season increases vegetation exposure to frost[END_REF][START_REF] Ma | Divergent trends in the risk of spring frost damage to trees in Europe with recent warming[END_REF][START_REF] Zohner | Late-spring frost risk between 1959 and 2017 decreased in North America but increased in Europe and Asia[END_REF] in spring, potentially resulting in dramatic impacts on agriculture [START_REF] Xiao | Estimating spring frost and its impact on yield across winter wheat in China[END_REF][START_REF] Unterberger | Spring frost risk for regional apple production under a warmer climate[END_REF][START_REF] Leolini | Late spring frost impacts on future grapevine distribution in Europe[END_REF] and forestry [START_REF] Greco | Late Spring Frost in Mediterranean Beech Forests: Extended Crown Dieback and Short-Term Effects on Moth Communities[END_REF][START_REF] Augspurger | warmth and frost: phenology, damage and refoliation in a temperate deciduous forest[END_REF] . The use of energy balances to study and better predict these environmental risks can provide novel insights into the responses of plants to extreme temperatures and offer more robust predictive tools, which are essential for developing adaptation measures and reducing the ecological and economic impacts. Temperature of plant organs and their dynamics are still overlooked in both environmental studies and modeling exercises [START_REF] Dong | Biophysical homoeostasis of leaf temperature: A neglected process for vegetation and land-surface modelling[END_REF] . Energy balance thus plays a key role, not only for plant phenology but also for all other processes since plant tissue temperature will govern key mechanisms such as photosynthesis and respiration and the general functioning of the plant.

Online methods

Description of the bud energy budget model

We implemented a simplified energy budget model for buds based on Landsberg et al. [START_REF] Landsberg | Apple bud and blossom temperatures[END_REF] , Hamer [START_REF] Hamer | The heat balance of apple buds and blossoms. Part I. Heat transfer in the outdoor environment[END_REF] , as well as equations from Jones 50 and Muir [START_REF] Muir | tealeaves: an R package for modelling leaf temperature using energy budgets[END_REF] . We refer to Supplementary Table S1 for parameter values and units. Our implementation follows the 'big-leaf' concept for which the whole canopy is simplified by a single representative organ, here an isolated bud in a deciduous canopy.

For an isolated bud, the amount of absorbed incoming radiation (Rabs) is balanced by the thermal infrared radiation loss (LWbud) and the energy lost by conduction and convection, generally called the sensible heat flux (H). The thermal time constant of buds is less than a minute [START_REF] Michaletz | A heat transfer model of crown scorch in forest fires[END_REF] , meaning that bud temperature responds quite rapidly to local changes. Because most phenology studies use meteorological data from local stations or gridded datasets, here we simulate bud temperature by considering that the energy balance is close to equilibrium at a time scale of a few minutes, which is consistent with temporal resolution of meteorological observations from FLUXNET sites (30 min) and CRU-JRA 68 (6 h) and the thermal time constant of buds ranging from a few seconds to a few minutes [START_REF] Michaletz | A heat transfer model of crown scorch in forest fires[END_REF] :

𝑅 𝑎𝑏𝑠 = 𝐿𝑊 𝑏𝑢𝑑 + 𝐻 ( 1 )
As we consider the energy budget of the whole, two-sided bud, and not the projected area, the amount of absorbed energy by buds is the sum of incoming shortwave (SW, visible and near-infrared) and longwave (LW, infrared) radiations from the sky and the ground:

𝑅 𝑎𝑏𝑠 = 𝛼 𝑠𝑤 (1 + 𝑟)𝑆𝑊 + 𝛼 𝑙𝑤 (𝐿𝑊 𝑠𝑘𝑦 + 𝐿𝑊 𝑔𝑛𝑑 ) ( 2 )
where αsw is the bud absorptivity to SW; r is the fraction of SW reflected by the ground and αlw is the bud absorptivity to LW, which is here defined as the average of LW coming from the atmosphere LWsky and coming from the ground LWgnd.

LW emitted by surrounding objects such as branches were considered equal to LW emitted from the ground and the sky. We set αlw to 0.97, which corresponds to the average absorptivity for wood and leaves [START_REF] Jones | Plants and microclimate. A quantitative approach to environmental plant physiology[END_REF][START_REF] Dupleix | Mid-infrared absorption properties of green wood[END_REF][START_REF] Howard | IR Radiation from Trees to a Ski Run: A Case Study[END_REF] . Since no data were available in the literature, we tested two different values of αsw, 0.5 and 0.8, corresponding to values commonly used for broadleaves and needleleaves, respectively. We set r to 0.2, which corresponds to a reasonable value for the fraction of SW reflected by grasses. Of course, r will strongly vary according to the albedo of the ground (e.g. understory/grass, forest litter, snow, etc.), which was simplified here for our perspective paper.

LWgnd was computed from ground temperature following the Stefan-Boltzmann equation:

𝐿𝑊 𝑔𝑛𝑑 = 𝜎𝜀𝑇 𝑔𝑛𝑑 4 ( 3
)
where σ is the Stefan-Boltzmann constant and ε is the bud emissivity to longwave radiations (which is equal to αlw since plant material tends to behave like a black body in the IR spectrum [START_REF] Jones | Plants and microclimate. A quantitative approach to environmental plant physiology[END_REF] ).

This is important to keep in mind that ground temperature will strongly depend on soil type, vegetation and snow cover but also soil humidity.

Buds lose thermal infrared radiation proportionally to their temperature as:

𝐿𝑊 𝑏𝑢𝑑 = 2𝜎𝜀𝑇 𝑏𝑢𝑑 4 ( 4
)
where σ is the Stefan-Boltzmann constant and ε is bud emissivity to longwave radiations. As for the absorbed radiation, radiations are emitted on both sides of the bud (i.e. toward the sky and the ground).

Finally, the sensible heat flux depends on the air to bud temperature gradient and is formulated as in Muir [START_REF] Muir | tealeaves: an R package for modelling leaf temperature using energy budgets[END_REF] :

𝐻 = 2𝜌 𝑎 𝑐 𝑝 𝑔 𝑏 (𝑇 𝑏𝑢𝑑 -𝑇 𝑎𝑖𝑟 ) ( 5 )
where ρa is the density of the air; cp is the specific heat capacity of air at constant pressure and gb is the boundary-layer conductance to heat.

𝜌 𝑎 = 𝑃 𝑅 𝑎 𝑇 𝑎 (1 + 𝑞 𝑎 ) (1 + 𝑅 𝑤 𝑅 𝑎 𝑞 𝑎 ) ( 6
)
where P is the atmospheric pressure, Ra and Rw are the specific gas constant for dry air and water vapor, respectively, Ta is the temperature of the air and qa is the specific humidity of the air.

𝑔 𝑏 = 𝐷 ℎ 𝑁𝑢 𝑑 ( 7
)
where Nu is the Nusselt number, Dh is the diffusion coefficient of heat in air and d is the bud diameter.

𝐷 ℎ = 𝐷 ℎ,0 ( 𝑇 273.15 ) 𝑒𝑇 101.3246 𝑃 ( 8 )
Dh is function of temperature and pressure, Dh,0 corresponds to Dh at 0°C and eT is the temperature dependence of diffusion.

The Nusselt number is estimated as a mixed convection such as:

𝑁𝑢 3.5 = 𝑁𝑢 𝑓𝑜𝑟𝑐𝑒𝑑 3.5 + 𝑁𝑢 𝑓𝑟𝑒𝑒 3.5 ( 9 )
with 𝑁𝑢 𝑓𝑜𝑟𝑐𝑒𝑑 = 𝑒 + 𝑎𝑅𝑒 𝑏 ( 10 ) and

𝑁𝑢 𝑓𝑟𝑒𝑒 = 𝑓 + 𝑐𝐺𝑟 𝑑 ( 11
)
Re and Gr are the Reynolds and Grashof numbers, respectively; a, b, c and d are constants that are dependent of the flow regime.

In our example, we simplified bud's shape as a spherical object with a diameter of 5 mm to computes the boundary-layer conductance to heat. Please refers to Hamer [START_REF] Hamer | The heat balance of apple buds and blossoms. Part I. Heat transfer in the outdoor environment[END_REF] for a detailed discussion of this assumption, as well as empirical formulations for the convective heat transfer of apple buds. Condition for laminar and turbulent flows, as well as constants a, b, c and d for a spherical object, as well as for objects of various forms, can be found in Monteith and Unsworth [START_REF] Monteith | Principles of environmental physics. Plants, animals, and the atmosphere[END_REF] . This is a simplification to keep in mind when calculating convection processes for different species with various bud shapes and for which the close proximity of branches can modify the flow regime. Also, the evolution of turbulence and sensible heat transfer from the bud to the canopy scale (associated with aerodynamic conductance in addition to bud boundary layer one) was not considered here, and was implicitly accounted for in the form of air temperature.

In addition, we performed two simulations varying ground albedo from 0.1 (~wet bare soil) to 0.9

(~snow) and bud diameter from 5 to 13 mm (Fig. 5 and6)

Note that a proper calibration of species and site properties (e.g. albedo) is needed to use this model for predictions. Here, we used already existing models and assumptions for the only purpose of exploring the expected variability in bud temperature with environmental conditions and for an isolated object representative of apical buds in a tree. Model development and validation was not intended here and will require in situ observations of bud traits and temperature data.

The bud energy model simulates bud temperature in idealized conditions and at equilibrium, without distinction between species and sites and several simplifications as described above.

Still, we describe in the two following sections how to account for bud energy storage as well as latent heat from wet buds in rainy conditions. Both implementations were tested without significant changes in the results and the corresponding code is also implemented and available in the model.

Accounting for energy storage

We applied the same approach accounting for the thermal time constant of buds, which is dependent on bud traits and its associated thermal resistance. Following Jones [START_REF] Jones | Plants and microclimate. A quantitative approach to environmental plant physiology[END_REF] , the change in bud temperature with time (t in s) follows:

𝑑𝑇 𝑑𝑡 = 𝑆 𝜌 * 𝑐 𝑝 * 𝑙 * ⁄ ( 12 )
with S the energy stored by the bud in W m -2 , ρ * and cp * are the density in kg m -3 and the specific heat capacity in J kg -1 K -1 of the bud and l * is the volume to area ratio, corresponding to d/4 for a cylinder and d/6 for a sphere, with d the bud diameter in m.

Between two equilibrium states, this first-order differential equation gives the following bud temperature:

𝑇 𝑏𝑢𝑑 (𝑡) = 𝑇 𝑒𝑞 (𝑡) -[𝑇 𝑒𝑞 (𝑡) -𝑇 𝑏𝑢𝑑 (𝑡 -1 )]exp (-𝑡 𝜏 ⁄ ) ( 13 )
with Teq the bud temperature at equilibrium and τ the thermal time constant in s (see eq. 9.8, eq. 9.9 and eq. 9.10, p. 227 in Jones [START_REF] Jones | Plants and microclimate. A quantitative approach to environmental plant physiology[END_REF] ) which is function of the conditions in the boundary layer [START_REF] Michaletz | A heat transfer model of crown scorch in forest fires[END_REF][START_REF] Bergman | Fundamentals of heat and mass transfer[END_REF] (i.e. bud traits and wind):

𝜏 = 𝜌 * 𝑐 𝑝 * 𝑙 * ℎ𝐴 ( 13 )
with A the wetted area in m² and h the heat transfer coefficient in W m -2 K -1 following:

ℎ = 𝑁𝑢𝜆 𝑎 𝑑 ( 13 )
with λa the thermal conductivity of air in W m -2 K -1 .

Michaletz & Johnson 55 estimated average values of τ ranging from 8 to 30 s for 12 different species with an average value of 13 s. The effect of energy storage is thus negligible with >30 min resolution meteorological data. Using high resolution micrometeorological observations should lead to better results by accounting for rapidly changing radiation load with clouds, wind speed, bud inclination etc.

Accounting for wet conditions

We implemented an option to account for latent heat in the model in order to simulate the cooling effect of evaporation when buds are wet because of the rain or dew conditions.

The maximum water density (dw) on buds depends on bud orientation, size and hydrophilicity. Since no data were found in the literature regarding water interception from buds, we defined the maximum interception reservoir based on leaf data. In the literature, dw for leaves varies between 0.05 and 0.2 mm 73 .

We thus arbitrary set dw for buds to 0.15 mm. Varying dw from 0 to 0.2 mm (kg m -2 ) did not change the results.

The maximum surface water content Wcmax (in kg) intercepted by a bud was thus calculated as:

𝑊𝑐 𝑚𝑎𝑥 = 𝑑 𝑤 × 𝑆 𝑏 ( 14 )
With Sb the surface area of the bud in m².

Surface water content (Wc) was estimated at each time step as the balance between water inputs from rain and dew and water outputs from evaporation. First, water coming from the rain of the previous time step is added to current Wc as:

𝑊𝑐 = 𝑚𝑖𝑛(𝑊𝑐 𝑚𝑎𝑥 , (𝑊𝑐 + 𝑟𝑎𝑖𝑛 × 𝑆 𝑏 )) ( 15 )
Then, latent heat energy dissipated by evaporation or accumulated from condensing dew (E) was calculated based on Wc following Gerlein-Safdi et al. [START_REF] Gerlein-Safdi | Dew deposition suppresses transpiration and carbon uptake in leaves[END_REF] :

𝐸 = { 𝜆 𝑣 (𝑇 𝑏𝑢𝑑 ) 𝑑𝑊𝑐 𝑑𝑡 , 𝑖𝑓 𝑑𝑊𝑐 𝑑𝑡 > 0 𝜆 𝑐 (𝑇 𝑏𝑢𝑑 ) 𝑑𝑊𝑐 𝑑𝑡 , 𝑖𝑓 𝑑𝑊𝑐 𝑑𝑡 < 0 ( 16 )
with the latent heat of condensation λc (in J kg -1 ) being equal of the latent heat of vaporization λv. The surface water content balance can be written as:

𝑑𝑊𝑐 𝑑𝑡 = 0.622 𝑆 𝑏 𝜌 𝑎 𝑔 ℎ [ 𝑒 𝑐 (𝑇 𝑎𝑖𝑟 ) -𝑒 𝑠𝑎𝑡 (𝑇 𝑏𝑢𝑑 ) 𝑃 ] ( 17 )
with gh the conductance of the boundary layer to water vapor in m s -1 , ec the vapor pressure in Pa and esat the saturating vapor pressure in Pa.

Here, water does not modify the diffusion coefficient of heat in air Dh.

Accounting for latent heat in the energy budget did not change the results.

FiguresFigure 1 |

 1 Figures

Figure 2 |

 2 Figure 2 | Simulated differences in temperature between buds and the air (ΔT) from energy balance. a. Daily variation in ΔT for an exposed bud and a typical day in April for two solar absorptivities: α=0.5 (typical for broad leaves 50 ) and α=0.8 (typical for needle leaves 50 ). SW and LWsky correspond to the amount of incoming shortwave and longwave radiation from the sky. b Example of ΔT simulated under idealized conditions (e.g. exposed buds in a deciduous canopy) using meteorological observations for winter and spring collected at the Hainich FLUXNET site 51 (temperate deciduous forest).The grey lines represent half-hourly differences in temperature simulated using an energy-budget model. The blue, black and red curves represent the 10-d rolling mean of the minimal (ΔTmin), average (ΔTave) and maximal (ΔTmax) temperature differences, respectively. Approximative Start Of growing Season (SOS) is illustrated by the light red region.

Figure 3 |

 3 Figure 3 | Spatial variability in bud temperature. Temperature difference (ΔT in °C) between buds and air simulated for four different FLUXNET sites during the year 1998: DK-Sor (Deciduous Broadleaf Forest; 55.48 °N, 11.64 °E), CA-Oas (DBF; 53.63° N, 106.20° W), US-Ha1 (DBF; 42.53 °N, 72.17 °W) and IT-Col (DBF; 41.85°N, 13.59°E). The grey lines represent half-hourly differences in temperature simulated using the energy-budget model. The blue, black and red curves represent the 10-d rolling mean of the minimal (ΔTmin), average (ΔTave) and maximal (ΔTmax) temperature differences, respectively. Approximative Start Of growing Season (SOS) based on fluxes is illustrated by the light red region. Note that these simulations represent idealized conditions and use the same parameterization without distinction between sites and species.

Figure 4 |

 4 Figure 4 | Potential changes in Tbud-Tair over Europe. a. Each point corresponds to the mean preseason ΔTave (see Fig. 2) simulated for six deciduous species across Europe (1059 sites) under idealized conditions (i.e. sun-exposed buds) using field observation of budburst and global meteorological data (see Supplementary Material). All sites and species were pooled together. Two solar absorptivity values were tested, 0.5 (grey) and 0.8 (black). The error bars represent the spatial and species variability (±1 SD around the mean, n = 5050). b. Under these conditions, ΔT is expected to decrease (blue) over the period 1990-2015 (n = 375 sites*species, 7%) and increase (n = 670 sites*species, 13% in orange). The error bars represent the spatial and species variability (±1 SD around the mean). The black points correspond to all sites*species pooled together (n = 5050) with α = 0.8 as illustrated in panel a.

Figure 5 |Figure 6 |

 56 Figure 5 | Impact of ground albedo on simulated bud temperature. Temperature difference (ΔT in °C) between buds andair simulated for two different ground albedo, r = 0.1 (~wet bare soil) and r = 0.9 (~snow). Temperature difference was simulated using meteorological observations for winter and spring collected at the Hainich FLUXNET site (temperate deciduous forest, Germany; see Figure2). The differences for maximum (red), average (black) and minimum (blue) temperatures are illustrated. Here, bud solar absorptivity to shortwave radiation and bud diameter were set to 0.8 and 5 mm, respectively

Acknowledgments

The authors would like to acknowledge Dr. Chris Muir and Dr. Renée Marchin Prokopavicius for their constructive feedbacks on our work. M.P. would like to acknowledge the financial support from the Fonds Wetenschappelijk Onderzoek (FWO; grant no. G018319N) and the H2020 Marie Skłodowska-Curie Actions (LEAF-2-TBM grant no. 891369). J.P. would like to acknowledge the financial support from the European Research Council Synergy grant ERC-SyG-2013-610028 IMBALANCE-P, the Spanish Government grant PID2019-110521GB-I00, the Fundación Ramon Areces grant ELEMENTAL-CLIMATE, and the Catalan Government grant SGR 2017-1005. H.V. acknowledges the support from the European Research Council Starting Grant 637643 TREECLIMBERS.

Data availability FLUXNET2015 data are available at https://fluxnet.fluxdata.org/data/fluxnet2015-dataset/ CRU-JRA data are available at https://catalogue.ceda.ac.uk/uuid/13f3635174794bb98cf8ac4b0ee8f4ed ERA5 soil temperature data are available at https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land PEP725 phenology data are available at http://www.pep725.eu/ Code availability The R code of the model of energy budgets and data sets used to generate the figures and analysis of this manuscript are available from Github at https://github.com/mpeaucelle/Tbud A version of the git repository is archived on Zenodo at https://zenodo.org/record/5897267 corresponding to tag v.2.0.

Datasets and analysis

Half-hourly forcing meteorological and soil temperature data (Fig. 2, 3, 5 and 6) were downloaded from the FLUXNET2015 dataset at https://fluxnet.fluxdata.org/data/fluxnet2015-dataset/. Historical climate data from CRU-JRA [START_REF] Harris | CRU JRA v1.1: A forcings dataset of gridded land surface blend of Climatic Research Unit (CRU) and Japanese reanalysis (JRA) data[END_REF] at a spatial resolution of 0.25° and at a temporal resolution of 6 h (Fig. 4 and Extended Fig. 2) were downloaded at https://catalogue.ceda.ac.uk/uuid/13f3635174794bb98cf8ac4b0ee8f4ed. Soil temperature data at 0.25°

were downloaded from ECMWF-ERA5 75 reanalysis with the 'KrigR' package [START_REF] Kusch | KrigR --A tool for downloading and statistically downscaling climate reanalysis data[END_REF] .

In situ leaf unfolding data (Fig. 4 and Extended Fig. 2) for Common alder (Alnus glutinosa), horse chestnut (Aesculus hippocastanum), silver birch (Betula pendula), European beech (Fagus sylvatica),

European ash (Fraxinus excelsior) and pedunculate oak (Quercus robur) were downloaded from the Pan European Phenology network (http://www.pep725.com/). Phenological observations followed the Biologische Bundesanstalt, Bundessortenamt und Chemische Industrie (BBCH) code, with leaf unfolding corresponding to BBCH = 11. Only sites with more than 20 years of observation over the 1990-2015 period were used, corresponding to 5050 sites*species in total, covering 1059 sites (Extended Fig. 1). Analysis and figures were generated with the R v3.5.1 software [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] .

Author Contributions Statement

M.P. designed the study, performed the analysis and wrote the first version of the manuscript. J.P. and H.V. substantially contributed to the interpretation of the results and the revisions of the manuscript. All authors read and approved the final manuscript.

Competing Interests Statement

The authors declare no competing interests


domain: sdv stat

text:
Sarah Berdot
email: sarah.berdot@aphp.fr

Cécile Blanc

Delphine Chevalier

Yvonnick Bezie

Laetitia Minh Maï Lê

Brigitte Sabatier

Impact of drug storage systems: a quasi-experimental study with and without an automated-drug dispensing cabinet

Keywords: healthcare costs, hospital care, information technology, medication process errors, pharmacy

Objective: To compare the costs and benefits of an automated-drug dispensing cabinet (ADC) versus traditional floor stock storage (TFSS). Design: A quasi-experimental multicenter study conducted during 2015. Setting: A teaching hospital (814 beds) equipped with 43 ADCs and a not-for-profit teaching hospital (643 beds) equipped with 38 TFSS systems, in Paris, France. Participants: All the wards of the two hospitals were included in the study. Intervention(s): ADC versus TFSS. Main Outcome Measure(s): A composite outcome composed of cost and benefits. Results: The total cost with payback period was substantially higher for the ADCs (574 006€ for 41 ADCs) than TFSS (190 305€ for 30 TFSS systems). The mean number of costly drugs and units were significantly higher for ADCs (P < 0.001). There was no significant difference in the mean number of overall drugs and units. There were significantly fewer urgent global deliveries with ADCs than TFSS units. Nurses' satisfaction with ADCs was high and the prevalence of medication process errors related to ADCs was low. No event due to storage errors was reported for ADCs and nine events were reported for TFSS units. On the contrary, informatic-related events increased with the use of ADCs, as expected. Conclusions: Overall, ADCs are well-established in wards and are particularly appreciated by nurses. A significant difference in the initial investment cost was confirmed, but it must be adjusted over time. This difference is offset in the long-term by gains in preparation time and fewer medication process errors, securing the medication process.

Introduction

Automated-drug dispensing cabinets (ADC) have been used since the 1980s in North America, but their introduction in France is more recent, dating to the late 1990s. In Canada, 36% of hospitals responding to a survey in 2010 were equipped with ADCs and this rate increased to 50% in 2014 [START_REF]The Hospital Pharmacy in Canada Editorial Board[END_REF]. Pedersen et al. estimated that 97% of responding hospitals in the USA were equipped with ADCs in 2014 [START_REF] Pedersen | ASHP national survey of pharmacy practice in hospital settings: Dispensing and administration-2014[END_REF], compared with 58% in 2002 [START_REF] Pedersen | ASHP national survey of pharmacy practice in hospital settings: dispensing and administration-2002[END_REF]. The floor stock management and rationalization induced by ADCs is important for North American countries where the system of unit billing to the patient is available. In France, the organization of drug dispensing in hospitals is still heterogeneous in terms of equipment (dispensing robot in pharmacy, ADCs, floor stock storage) and stock management (global and unit-dose dispensing, with urgent global deliveries for specific needs). A report by the French hospitalization and care organization directorate (DGOS) in 2009 found that only 2% of French hospitals (any type of healthcare institute) had deployed one or more ADCs [START_REF]Etude d'impact organisationnel et économique de la sécurisation du médicament dans les établissements de santé[END_REF]. The total number of ADCs in French healthcare institutes represented 0.3% of the total number of pharmacy storage systems. This low rate can be explained by the high unit cost of ADCs, without evaluation of the return on investment and limited manufacturers and products available in France.

Numerous studies have been carried out to measure the impact of ADCs on the medication process in clinical wards [START_REF] Borel | Effect of an automated, nursing unit-based drugdispensing device on medication errors[END_REF][START_REF] Cousein | Effect of automated drug distribution systems on medication error rates in a short-stay geriatric unit[END_REF][START_REF] Chapuis | Automated drug dispensing system reduces medication errors in an intensive care setting[END_REF][START_REF] Risør | Complex automated medication systems reduce medication administration errors in a Danish acute medical unit[END_REF]. Most evaluated the impact following implementation in a single hospital. They showed a potential positive impact on the medication process by reducing medication errors at every step of the process, from filling, picking and preparation to administration. Ray et al. reported a 32% decrease in refill errors [START_REF] Ray | Experience with an automated point-of-use unit-dose drug distribution system[END_REF]. Administration errors also decreased from 53% to 39% with the use of ADCs [START_REF] Borel | Effect of an automated, nursing unit-based drugdispensing device on medication errors[END_REF][START_REF] Cousein | Effect of automated drug distribution systems on medication error rates in a short-stay geriatric unit[END_REF]. A randomized, controlled study showed a significant decrease in the prevalence of medication errors (18.6% with traditional floor stock storage (TFSS) and 13.5% with ADCs, P < 0.05) [START_REF] Chapuis | Automated drug dispensing system reduces medication errors in an intensive care setting[END_REF]. However, studies focusing on the economic impact of ADCs are few or outdated [START_REF] Chapuis | Automated drug dispensing systems in the intensive care unit: a financial analysis[END_REF][START_REF] Lee | Use of an automated medication storage and distribution system[END_REF][START_REF] Schwarz | Implementation and evaluation of an automated dispensing system[END_REF]. In a recent economic analysis in France, Chapuis et al. [START_REF] Chapuis | Automated drug dispensing systems in the intensive care unit: a financial analysis[END_REF] found that the implementation of ADCs was financially profitable, with a positive payback period of 5 years.

The European Georges-Pompidou Hospital has progressively invested in ADCs since its opening in 2000 to optimize and secure the medication process. The benefits of these cabinets in the pharmacy and clinical wards merit evaluation, given their cost. The objective of this study was to compare the cost and benefits of dispensing systems in two hospital settings, one with ADCs and the other with TFSS.

Methods

Population and setting

This study was conducted in two French hospitals. The first, the European Georges-Pompidou Hospital (HEGP), is an 814-bed teaching hospital in Paris, France. HEGP has 43 ADCs (OmniRx; Omnicell, Mountain View, CA, USA). The ADCs were implemented from 2000 to 2014. In each ward, all drugs are stored in the ADCs, and nowhere else. Drugs are stored in a random order. The cabinets are locked, and the nurses have access through passwords or fingerprint identification. They have access to all drugs, narcotics included. After selecting the patient name, nurses select the drug prescribed from the global list. The prescription software is not yet interfaced with the ADC software. The ADC software records the time, the identity of the nurse, the patient's name and the drugs removed. The ADC compartments are refilled weekly by pharmacy technicians and the quantities of drugs determined using the interfaced pharmacy information system to assess the need.

The second hospital is the Paris Saint-Joseph hospital (HPSJ), a 643-bed not-for-profit teaching hospital located in Paris, France. Drugs are stored in 38 TFSS systems from Herman Miller ® , which have been in use since 1990. These open floor stock systems consist of one or two cabinets. Storage is in alphabetical order, with a separation between oral and injectable drugs. Narcotic drugs are managed outside the system in locked cabinets. The TFSS units are refilled weekly by pharmacy technicians, as for the ADCs.

Both hospitals are equipped with a patient information system integrating an electronic health record and a computerized physician order entry system (DxCare ® , Medasys). The pharmacy software used is Pharma ® (Computer Engineering). Global deliveries using Pharma ® are scheduled for weekly refilling in both hospitals. This refilling normally satisfies urgent needs. However, it is sometime necessary to complete the stock with punctual urgent global deliveries.

Study design

We conducted a quasi-experimental multicenter study from January to December 2015.

Inclusion criteria

All the wards of the two hospitals were included in the study, including those with an average length of stay above 24 h, emergency departments and intensive care units.

Outcomes

The main outcome was a composite outcome composed of cost and benefits calculated for the study period (Table 1). The cost criteria were derived from Chapuis et al. [START_REF] Chapuis | Automated drug dispensing systems in the intensive care unit: a financial analysis[END_REF] and consisted of (1) the cost of storage, (2) the cost of the immobilized drug stock and (3) the cost to employer of pharmacy technicians. The benefits criteria consisted of (1) the number and quantity of each drug in the storage, (2) the number of refillings and drug lines reloaded, (3) the cumulative refiling time, (4) nurses' satisfaction and (5) medication process errors declared during the study period.

Statistical analysis

Data are expressed as the mean and standard deviation (SD) for quantitative variables and as the number and associated percentage for qualitative variables. Group comparisons were performed using the Wilcoxon test for quantitative variables and the Student test for qualitative variables. Data were analyzed with R software version 3.3.0. Satisfaction was compared with the Chi-square test and exact Fisher test using BioStatgv (http://marne.u707.jussieu.fr/biostatgv/).

Results

For this study, 23 wards in HEGP (41 ADCs) and 21 in HPSJ (30 TFSS systems) were included. The total number of hospitalizations was similar for both hospitals (33 351 for HEGP vs. 37 070 for HPSJ). The average length of stay was 10.2 days (min: 2.7-max: 18, mean of 30 beds/ward) in HEGP and 4.6 days (min: 1.1-max: 8.9, mean of 26 beds/ward) in HPSJ. The mean age of the patients was 61 years for both hospitals.

Cost

We evaluated the costs during the study period, integrating the initial costs of investment with the payback period and those following the payback period (Table 2). The mean cost with and following the payback period for one storage was significantly different between the two systems (P < 0.001).For the ADCs, the investment costs were 21 840€ for simple ADCs and 24 340€ for double ADCs, with software, license and annual maintenance costs. The cost of immobilized drug stock was 194 805€ for 2015. For TFSS, the investment costs were 3500€ and 7000€ for one and two cabinets, respectively. There were no associated software or license costs. The maintenance costs could not be evaluated. The cost of immobilized drug stock For both storage systems (ADCs or TFSS), the manufacturers estimated an amortization period of 5 years, with a repartition of the total investment costs (storage, software, and licenses between the pharmacy department and the ADC computer system) over this period. Thus, during the first 5 years, the cost of storage took into account 20% of the investment costs and the cost of annual maintenance (as storage deployment was continuous). We also evaluated the case after this payback period, in which only annual maintenance costs were considered to calculate the storage costs. (2) Cost of immobilized drug stocks Cost of immobilized drug stocks for 1 year (2015)

The cost of immobilized drug stocks was estimated by floor stock inventories. Nurses' satisfaction The nurses' satisfaction was evaluated through a self-administered questionnaire derived from a validated questionnaire developed by Palm et al. [START_REF] Palm | Determinants of user satisfaction with a Clinical Information System[END_REF][START_REF] Palm | Determinants of clinical information system post-adoption success[END_REF]. Four dimensions were estimated (compatibility of the system, usefulness, ease of use and user satisfaction) for the ADC and TFSS systems. A three-point Likert scale was used: satisfied, unsatisfied and not satisfied or unsatisfied (5) Medication process errors declared in 2015

Medication process errors Medication process errors were analyzed using the medication process errors registry available in each hospital, in which every healthcare professional can declare any medication process errors. All events related to the medication process declared during the studied period were evaluated. was 52 405€. The mean cost of immobilized drug stock for one storage was significantly lower than the same cost for ADCs (P < 0.001). Finally, the total cost, with payback period, was 574 006€ for ADCs and 190,305€ for TFSS systems. Following the payback period, the total cost for one storage system was 9237€ for ADCs and 5106€ for TFSS.

Benefits

There was no significant difference in the mean number of drugs (P = 0.79) or drug units (P = 0.07) stored in the systems (Table 3). The mean number of costly drugs and units were significantly higher for ADCs (P < 0.001). The mean number of global deliveries for ADCs was significantly higher than that for TFSS (99 vs. 72, P < 0.001). There were significantly fewer urgent global deliveries with the ADCs than with TFSS. The number of urgent global deliveries represented 70% of the total requests (global deliveries and urgent global deliveries) for ADCs and 90% for TFSS. Thirty self-measurements were collected at HEGP and 13 at HPSJ. The average time per refilling line obtained was 1.65 min for ADCs and 1.60 min for TFSS. Based on the 2 days of observation performed by the observer (CB), the times were 1.66 min for ADCs and 1.67 min for TFSS. There was no difference between the two systems. The total annual times for dispensing are presented in Table 3.

Fourteen of 23 (61%) wards answered the user satisfaction questionnaire at HEGP and 14 of 21 (67%) at HPSJ. Fifty-six questionnaires from HEGP and 25 from HPSJ were collected. Four nurses in HPSJ had experience with both storage systems and 75% preferred ADCs. In some wards, several nurses and nurse managers participated. Satisfaction rates are presented in Table 4. Overall, user satisfaction was significantly higher with ADCs than TFSS (89% vs. 56%, P < 0.001). Despite the perception of comparative compatibility of both storage systems (P = 0.06), there were significant differences for the ease of use (P < 0.001) and usefulness dimensions (P < 0.001), with better perception of the ADC system.

During the study period, 75 medication process errors related to a 'drug event: medication error process/dysfunction' were reported at HEGP. Among these, 12 were related to ADCs (16%). At HPSJ, there were 40 errors related to TFSS among 99 events reported (40%). Four types of medication process errors were identified: storage refill management, informatics problems, administration errors and other types of events (Fig. 1). Storage refill management events were reported five times more often for TFSS than for ADCs. There were seven events reported at HEGP (six due to treatment supply failure and one due to an expired drug present in the ADC). At HPSJ, there were 19 events due to treatment supply failure, one to the presence of an expired drug, four linked to the lack of drug storage by pharmacy technicians, nine to storage of the wrong drug and six to delivery of the wrong drug. There were four events related to informatics problems, such as the lack of a nurse password or an error in the password to enter the ADC. Finally, there was one event reported at HEGP which generated a wrong administration error due to confusing one drug with another. The other event at HPSJ concerned a problem with refrigerator temperature.

Discussion

The total cost, with payback period, was 574,006€ for 41 ADCs and 190,305€ for 30 TFSS systems. The mean number of costly drugs and units were significantly higher for ADCs (P < 0.001). However, there was no significant difference in the mean number of stored drugs or units. There were significantly fewer urgent global deliveries with ADCs than TFSS. The total annual time for deliveries was 115 h for ADCs and 164 h for TFSS. Nurses' satisfaction with ADCs was high and the prevalence of medication process errors related to ADCs lower than that for TFSS. No medication process error due to storage errors was reported for ADCs, whereas nine were reported for TFSS. On the contrary, informatics-related events increased with the use of ADCs, as expected. Some of our results are consistent with those found in the literature [START_REF] Tsao | Decentralized automated dispensing devices: systematic review of clinical and economic impacts in hospitals[END_REF]. The use of ADCs at Robert Ballanger Hospital Center (France) resulted in fewer urgent global deliveries with a 75% decrease of demands during daytime care and 50% during night-time care [START_REF]Chedhomme FX Apport des armoires robotisées dans la sécurisation et l'optimisation du circuit du médicament[END_REF]. Chapuis et al. reported a 43-57% decrease in treatment failures in the two wards equipped with ADCs [START_REF] Chapuis | Impact d'une armoire à pharmacie sécurisée en réanimation médicale sur l'iatrogénèse médicamenteuse, la[END_REF]. User satisfaction with the use of ADCs has been evaluated in several studies [START_REF] Ottino | Impact des cabinets automatisés décentralisés sur la diminution des erreurs de médication en soins intensifs[END_REF]. In France, 78% of nurses have been reported to be satisfied or very satisfied with the daily use of ADCs [START_REF]Chedhomme FX Apport des armoires robotisées dans la sécurisation et l'optimisation du circuit du médicament[END_REF]. They appreciated the fact that the ADCs indicated the equivalents between brand name drugs and their generic term on the ADC screen and redirected them to the correct reference available in the system.

Our findings on reported medication process errors related to the storage system are consistent with the literature. Studies showed a decrease in the number of drug errors with ADCs, both in the refilling [START_REF] Ray | Experience with an automated point-of-use unit-dose drug distribution system[END_REF] and administration steps [START_REF] Borel | Effect of an automated, nursing unit-based drugdispensing device on medication errors[END_REF][START_REF] Cousein | Effect of automated drug distribution systems on medication error rates in a short-stay geriatric unit[END_REF]. In addition, three times fewer events due to treatment supply failures have been reported with ADCs. This observation is consistent with the higher number of urgent deliveries in HPSJ due to treatment failure.

In contrast, some of our results differ from those of the literature. Savings in the cost of immobilized drug stocks have been observed after implementing ADCs. These savings ranged from 29 to 56% [START_REF] Chapuis | Automated drug dispensing systems in the intensive care unit: a financial analysis[END_REF][START_REF] Sabatier | La mise en place d'un automate décentralisé dans une unité de réanimation médicale[END_REF][START_REF] Kheniene | Economic impact of an automated dispensing system in an intensive care unit[END_REF]. At HEGP, costly drugs were stored in the ADC, probably explaining this difference. The traceability and accuracy afforded by informatics has allowed the drugs stored in ADCs to be tracked, facilitating the storage of costly drugs and blood-derived drugs requiring patient traceability. Indeed, narcotics and blood-derived drugs were only stored in ADCs. Unlike in our study, several studies have reported an increase in the time of pharmacy technicians dedicated to ADC management [START_REF] Chapuis | Automated drug dispensing systems in the intensive care unit: a financial analysis[END_REF][START_REF] Tsao | Decentralized automated dispensing devices: systematic review of clinical and economic impacts in hospitals[END_REF][START_REF] Kheniene | Economic impact of an automated dispensing system in an intensive care unit[END_REF]. Some functions of ADCs allow better management and optimization of stocks. Urgent global deliveries were three times higher with TFSS when refilling drug stocks. This shows that the refilling planning is not adapted to a ward's needs and results in functional disorganization, even when a similar process is used by the two hospitals. The automation of the scanning step during the refilling process avoided a visit to the ward by pharmacy technicians and removed the need of transcriptions. Hospital stock visibility is especially appreciated during supply shortages, and the ADC itself avoided drug overstorage due to finite physical space allocated for a given drug. The risk of informatics failure should not be ignored. An override mode must be made available. Team training is essential for optimal use of ADCs, but it is more elaborate than that for TFSS. This is the first study to compare traditional storage and ADCs with a quasi-experimental design. Most studies have used a 'beforeafter' design to evaluate the implementation of ADCs within one or more wards. Our study has some limitations. The time-saving for nurses was not evaluated, but was predictable. Studies show a reduction in nurses' time spent on acquisition of a single narcotic dose and drug inventory, allowing more time spent with patients [START_REF] Tsao | Decentralized automated dispensing devices: systematic review of clinical and economic impacts in hospitals[END_REF]. There were two times fewer pharmacy technician selfmeasurements returned from HPSJ than HEGP. This difference can be explained by the less frequent presence of the pharmacist at HPSJ. More self-measurements and, more importantly, more observations would have made the results for refilling time even more robust. The medication process error reporting system is based on declarative reports. Underestimation of the prevalence of events is known for both storage systems. It would have been informative to directly observe users to evaluate medication errors using the two systems. Finally, it would have been informative to assess the loss of drug due to professional and personal consumption.

Overall, ADCs are well-established in wards and particularly appreciated by nurses. A significant difference in the initial investment cost was confirmed, which is offset in the long-term by gains in preparation time and fewer medication process errors, securing the medication process. The time of pharmacy technicians and nurses can be saved by minimizing urgent global deliveries and controlling their frequency, and functional disorganization may be avoided. The time saved can be used for ADCs maintenance by pharmacy technicians and patient care by nurses.

( 3 )

 3 Cost of pharmacy technicians Cost of pharmacy technicians for 1 year (2015) The cost of pharmacy technicians on medication-related work activities was calculated, considering an annual cost to employer of 42 000€ per technician. Benefits criteria (1) Number and quantity of each drug in storage Number of drugs and unit doses Number of costly drugs and unit doses For each system, the mean number of drugs and unit doses for costly drugs (more than 20€ per unit) and all drugs (costly and not costly) were calculated. (2) Number of refillings and drugs lines reloaded Number of global deliveries and urgent global deliveries Mean number of drug lines per global delivery The number of refillings was calculated by summing the number of global deliveries and urgent global deliveries. The mean number of drug lines per global delivery was also calculated based on all deliveries. (3) Cumulative time for refilling Cumulative time for refillingThe cumulative time for refilling was calculated by multiplying the total number of reloaded lines and the mean time for one reloaded line. The time for refilling one reloaded line was determined simultaneously by a pharmacy technician (self-measurement) and a pharmacist (CB). The average of the two times was used. All stages of refilling (scanning, management, preparation, storage and walks to the ward) were timed to obtain an average time per line. (4) Nurses' satisfaction

Table 1

 1 Description and measurement of the cost and benefits criteria

	Criteria	Description/measurement
	Cost criteria
	(1) Cost of storage	Included investment costs, the cost of purchasing software and licenses and annual maintenance costs

Table 2

 2 Costs for 41 ADCs and 30 TFSS systems

	ADCs (N = 41)	TFSS (N = 30)

*NE: Not evaluated.

Table 3

 3 Number, refilling and time of refilling for one storage system

	Type of storage	ADCs	TFSS	P-value
		Mean (SD)	Mean (SD)
	Number and quantity
	Mean number of drugs	228 (86)	233 (64)	0.79
	Mean number of costly drugs	3.5 (3.1)	1.4 (1.8)	<0.001
	Refilling
	Mean number of global deliveries (refilling)	99 (38)	72 (23)	<0.001
	Mean number of lines in global deliveries	3174 (1576)	4038 (1930)	<0.001
	Mean number of urgent global deliveries	226 (88)	677 (364)	<0.001
	Time of refilling
	Cumulative time for global deliveries	96 h	108 h
	Cumulative time for urgent global deliveries	19 h	56 h
	Total annual time for deliveries	115 h	164 h

Table 4

 4 Perception of drug storage (ADCs or TFSS) for the four dimensions of the user questionnaire-percent of satisfaction

		ADCs	TFSS	P-value
		N = 56	N = 25
	Compatibility*	93%	76%	0.06
	Usefulness*	98%	64%	<0.001
	Ease of use*	98%	66%	<0.001
	User satisfaction**	89%	56%	<0.001
	Statistical analyses were performed using the exact Fisher test (*) or Chi-
	square test (**).

Downloaded from https://academic.oup.com/intqhc/advance-article-abstract/doi/10.1093/intqhc/mzy155/5055358 by guest on 18 July 2018

Acknowledgements

The pharmacy technicians, pharmacy technician managers and nurse managers of both hospitals.

Funding

This study received no external funding.