Interannual Variations in Spring Snowmelt Timing of Alaskan Black Spruce Forests Using a Bulk‐Surface Energy Balance Approach

Spring snowmelt occurs for a short duration on an annual time scale, but their timings considerably affect the carbon and hydrological cycle in high‐latitude ecosystems. Here, we developed a simple snowmelt model, treating the ecosystem surface as a bulk‐surface layer. Energy fluxes across this bulk surface and the snow‐soil boundary determine snow temperature and the energy utilized for snowmelt. Parameterizing the bulk surface using decade‐long eddy covariance site data from two Alaskan open black spruce forests offered an opportunity to quantitatively evaluate meteorological drivers affecting snowmelt timings without the needs for detailed canopy information. The sensitivity analysis suggested that the total snowfall on the forest floor, ranging from 0.35 m in 2016 to about 1 m in 2018 and 2020, was the most crucial driver for snowmelt timing. This factor accounted for a 10‐day difference in the interannual variations in snow disappearance dates. The importance of the snowfall varied from year to year, and in 2013, the late snowmelt was characterized by low air temperatures, which increased sensible heat loss from the snowpack. The importance of atmospheric radiation was revealed in relatively warm years, such as 2016 and 2019. Our modeling approach necessitates adjusting one empirical parameter that reflects the heat conductivity from the bulk surface to the snow, based on observations. Nevertheless, despite this need for adjustment, the bulk‐surface approach helps identify important meteorological drivers underlying observed snowmelt within a simple theoretical framework.


Introduction
Spring snowmelt is an important hydrological event that affects the hydrological cycle and subsequent biological activity in high-latitude ecosystems.Spring snowmelt in interior Alaska typically commences in mid-April, and the entire landscape is snow-free within a few weeks.The timing of snow disappearance affects the onset of biological activity (Kawashima et al., 2021;Kobayashi et al., 2016;Suzuki et al., 2011;Zheng et al., 2022), and snowmelt water supports the evaporation demands of terrestrial biospheres during the growing season.Lateral water flows in the permafrost region provide carbon-rich water to riverine and coastal aquatic ecosystems (Ikawa & Oechel, 2014;Koch et al., 2014;Lougheed et al., 2020).Accurately evaluating the timing of snowmelts is therefore crucial for understanding ecosystem phenology, water resource availability, and terrestrial-aquatic ecosystem interactions in the high latitudes.
Snowmelt is a complex process that is influenced by various meteorological drivers, including temperature, radiation, wind speed, and humidity.Among these drivers, atmospheric radiation (longwave radiation flux) is considered a critical driver affecting snowmelt in high latitudes (Ohmura, 2001;Zhang et al., 1997), and the amount of its energy is generally greater than turbulent heat transfers (Hayashi et al., 2005;Ohmura, 2001;Suzuki et al., 2006).Humidity in the atmosphere can affect snowmelt through various mechanisms, including changes in atmospheric radiation (e.g., Iijima et al., 2007;Zhang et al., 1997) and latent heat transfer.Snowmelt timing is further complicated when the amount and distribution of snowfall are altered in the future (Hinzman et al., 2005;Serreze et al., 2002).Moreover, the key meteorological drivers that affect snowmelt would vary depending on the region and geographical characteristics (Kawashima et al., 2022;Langer et al., 2011;Mioduszewski et al., 2014;Sicart et al., 2008).Considering energy flows according to the energy balance theory would be helpful to disentangle individual effects on snowmelt.
While snow process models have been successful in evaluating snowmelt over flat snow surfaces (e.g., Bartelt & Lehning, 2002;Essery, 2015;Kominami et al., 2015;Niwano et al., 2012), the complexity of energy balance models increases when evaluating snow surface energy balance under plant canopies (Essery et al., 2009;Pomeroy et al., 2009;Rutter et al., 2009;Yamazaki et al., 2004).Radiation and thermal conditions inside the canopy structures of boreal forests are highly heterogeneous (Kobayashi et al., 2014;Nakai et al., 2023), making it further difficult to accurately describe the energy flows that determine spring snowmelt.Furthermore, Essery et al. (2009) highlighted the difficulty of parameterizations in snow process models.However, if our focus is solely on comprehending meteorological drivers underlying observed snowmelt and integrated ecosystem fluxes, rather than the processes within the canopy, a much simpler modeling approach may be helpful.Such an approach would take advantage of an ever-increasing amount of eddy covariance site data, which enables us to examine interannual variations of snowmelt without detailed canopy parameters.
This study aims to investigate the interannual variations in the spring snowmelt timings in open black spruce forests, which are the dominant ecosystem in interior Alaska.To this end, we developed a snowmelt model considering the ecosystem surface as a single bulk-surface layer.Assuming that the rate of snowmelt is determined by integrated energy fluxes occurring across the bulk surface and the snow-soil boundary, parameterizing the bulk-surface processes based on eddy covariance site data allows us to identify critical processes underlying observed snowmelt.The model was applied to the decadal data obtained from black spruce forest flux sites to identify the key meteorological drivers of the interannual variation in snowmelt timing.Additionally, we evaluated how different canopy surface properties, such as heat conductivity, roughness lengths, and albedo, affect snowmelt by comparing the canopy surface parameterizations of two open black spruce forests.This evaluation helps determine how varying canopy surface properties may impact snowmelt, even among the same biomes.

Black Spruce Study Sites
Two Alaskan black spruce flux sites were targeted for this study for the period from 2011 to 2020.In this study, we focused on the data collected in March-May  to analyze characteristics of spring snowmelt.The collaborative study site of the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) and International Arctic Research Center, University of Alaska Fairbanks (65°07′N, 147°29′W) is registered as US-Prr in the AmeriFlux network (Iwahana et al., 2023).Another black spruce forest site is located in the campus of University of Alaska Fairbanks (US-Uaf, 64°52′N, 147°51′W) (Ueyama et al., 2023).Black spruce (Picea mariana) is the dominant species over permafrost regions in boreal North America, and these two sites represent such a black spruce forest ecosystem in interior Alaska.The estimated leaf area index of the overstory during the growing season and canopy height are 0.73 m 2 m 2 (Kobayashi et al., 2014) and 2.91 m (Nakai et al., 2013), respectively, at the US-Prr site.The canopy structure of US-Uaf is very similar to US-Prr, with leaf area index that includes understory deciduous vegetation, ranges from 0.2 to 1.2 m 2 m 2 during the snowmelt season, and the canopy height is 3 m (Iwata et al., 2012).Snow density and conductivity surveys were also conducted at US-Uaf in 2011-2014 (Appendix A).At both sites, eddy covariance measurements have been operated for carbon fluxes and turbulent heat fluxes (i.e., sensible heat flux, H, and latent heat flux, λE) together with other meteorological observations.All data were half-hourly averaged.The reported energy imbalance in spring was 2% in US-Prr (Nakai et al., 2013) and 10% in US-Uaf (Iwata et al., 2012), and we did not apply any correction to H or λE related to the energy balance closure.
Table 1 describes a list of instruments used in this study.The eddy covariance system at US-Prr is positioned 11 m above the ground level and utilizes an enclosed infrared gas analyzer and a sonic anemometer.Snow depth was measured using a sonic ranging sensor at three locations, and the values were averaged.Soil temperature was measured using thermistors at six levels below a moss layer of about 0.14 m.Radiation fluxes were observed using a four-component radiometer with up-facing and down-facing pyranometers and pyranometers.Air temperature and relative humidity were measured at the same height as the eddy covariance system using a humidity and temperature probe installed in a ventilated shield.Further information is found in our previous works (Kobayashi et al., 2018(Kobayashi et al., , 2023;;Nagano et al., 2018).
At US-Uaf, the eddy covariance system was positioned 6 m above the ground, using an open-path infrared gas analyzer and a sonic anemometer.Radiation fluxes were observed using a four-component radiometer.Air temperature and relative humidity were measured with a humidity and temperature probe at the height of 8 m above the ground, enclosed in a ventilated shield.Continuous snow depth data from the winter of 2016 are available.For further information, the reader is referred to Ueyama et al. (2014) and Iwata et al. (2012).
For both sites, albedo (a) was calculated based on the daily average values of both downward and upward shortwave fluxes.Snowfall rate on the forest floor (P snow ) was estimated by calculating the temporal difference in the 24-hr running average of the snow depth data (d snow ).While rain-on-snow events alter d snow (McCabe et al., 2007), we considered such effects negligible at our site due to the consistently low air temperature when snow was present.

Definition of Snow Disappearance Date
We defined the observed snow disappearance date (SDD) as the day when the daily albedo became less than a certain threshold, because of the full data coverage of albedo data during the study period for both sites.The threshold albedo value of 0.3 was chosen using available snow depth data as the timing when daily average snow depth values fell below the accuracy of the sonic ranging sensor, that is, 10 mm (see Text S1 in Supporting Information S1).We also confirmed the determined SDDs corresponded to the dates when about half of the surface snow had disappeared, based on the phenology camera data from the US-Prr site between 2011 and 2016 (Nagai et al., 2013(Nagai et al., , 2018)).After this date, the surface snow disappeared in a few days.In the snowmelt model described in the next section, SDD was defined as the date when the snow water equivalent (d swe ) became less than 1 mm for the first time in spring.

Description
Figure 1 describes the framework of the model, and Table 2 lists the variables and constants used for the model.The time evolution of d swe was determined by the balance of the precipitation of snow (P swe ), the snowmelt (M snow ), and the rate of sublimation (E s ) (Equation 1).Then the balance may be written as follows: where ρ is the density, with its subscript denoting the property (air, snow, soil, or water).P swe was calculated from observed P snow and snow density (ρ snow ).E s was calculated from the modeled latent heat flux (λE) when M snow equaled zero.Otherwise, if snowmelt occurred (i.e., M snow > 0), E s was assumed to be zero.The depth of the snow d snow was calculated as One uniqueness of this model is that we introduce the ecosystem surface temperature (T eco ).T eco represents the temperature of the ecosystem surface, which is responsible for sensible heat exchange with the atmosphere as well as the upward longwave flux from the ecosystem surface to the atmosphere.The ecosystem surface during the targeted period of this study is a mixture of black spruce tree canopies and canopy gaps covered with snow.In this snowmelt model, we considered the snow as one layer regardless of the depth and assigned T snow as the representative temperature at the center of this layer.Furthermore, we regarded the downward heat flux from the ecosystem surface to the center of the snow layer (0.5d snow below the ecosystem surface) as G s .For the conductivity that determines G s , we introduced a term named bulk ecosystem surface conductivity (k eco ), which explains the complex effects of aboveground trees and snow surface interactions that occur inside the canopy.Using these variables, G s is written as  The temporal evolution of T snow is described as where c is the specific heat, and the subscript denotes the property (air, snow, or soil).The G 0 is the downward heat flux from the snow to the soil.The soil depth and temperature of the soil are denoted as z i and T soil,i .The i subscript is used when referring to a particular soil layer (i = 1, 2, … n where n = 7).The moss layer (0.14 m) was regarded as the first soil layer (i = 1), below which was grided at the depths of 0.0, 0.1, 0.3, 0.5, 0.7, 0.9, and 1.1 m.The grids were determined for the convenience of comparing with soil temperature data for US-Prr.The heat conductivity of the uppermost layer (k 1 ) was determined from the weighted average of the reciprocal values of heat conductivities of soil (k soil ) and snow (k snow ) along the distance between the centers of the snow layer and the uppermost soil layer (∆z 1 ), following the method used in Essery (2015).The λ f is the latent heat of fusion.The time evolution of ρ snow was estimated using a densification equation, using parameters estimated based on the snow-survey data in US-Uaf (Appendix A).The value of k snow was empirically determined from ρ snow .
The T soil of each layer was calculated using the following equations: where θ v is the porosity of the soil that was assumed to be fully saturated with water, and the value was fixed at 0.58, referring to volumetric water content data.ΔT f (≡ 1 K) is the temperature range in which the phase change occurs.The reasonability of the soil parameters (c soil ρ soil and k soil ) listed in Table 2 was checked with the time series of T soil at different depths (Heusinkveld et al., 2004) (see the detail in Text S2 in Supporting Information S1).The heat capacity and conductivity of the saturated and frozen moss layer were assumed to be the same as that of the soil.The top boundary condition of the soil layer was given by G 0 .
The bottom boundary condition T soil,n was estimated every 24 hours from the daily value of T soil,n 1 , assuming that their annual variations follow a sinusoidal curve (Hirota et al., 1995(Hirota et al., , 2002) ) as where D y is the damping depth, ω y [ = 2π/(365*86,400) s 1 ] is the angular frequency of the annual variation, and T m is the annual average of T soil,n (≡ 1°C).
The surface energy balance equations were solved for T eco and sensible and latent heat fluxes (H and λE, respectively) as where R sd is solar radiation (downward shortwave flux), R ld is atmospheric radiation (downward longwave flux), ε is the emissivity and absorptivity of the ecosystem surface (≡ 0.95), σ is Stefan-Boltzmann constant, U is wind speed, T air is air temperature, q sat (T eco ) is saturated specific humidity at T eco , and q air is specific humidity in the air.The bulk transfer coefficient for heat flux (C H ) was determined as the function of roughness lengths for the momentum and heat (z 0 and z T , respectively) as where κ is the von Kármán constant, and Ψ M and Ψ H are stability correction factors (Dyer & Hicks, 1970;Kondo et al., 1978) The observation height, denoted as z obs , has subscripts U or T indicating wind speed and temperature measurements, respectively.The zero-plane displacement d is fixed at 1.46 m according to the analysis in Chu et al. (2018).The evaporation efficiency (β) represents the moisture availability (Kondo et al., 1990) and corrects for the gradient of the specific humidity between the bulk surface and the air.We assumed that λE originated from sublimation (E s ) when M snow equaled zero.Otherwise, if snowmelt occurred, λE was considered to originate from evaporation.
The model was run from DOY 60 to SDD with initial values of observed d snow and T soil .The T soil data on DOY 60 in 2011 were not available, and their values were regarded as the same values in 2012.The model algorithm was coded with the Julia Programming Language version 1.9.2 (Bezanson et al., 2017).Figure 1b explains the computation flow as follows: 1. Compute the ecosystem surface energy balance (Equations 3 and 7) using the T snow value of the previous time step (half-hour) assuming the T snow change over 30 min was minimal.2. Solve T snow and T soil using the Crank-Nicolson scheme with the weighting factor of 0.6 (Equations 4 and 5), considering G s and T soil,n as the top and bottom boundary conditions, respectively.Equation 4 was first solved assuming M snow was zero.Then λ f M snow was computed as the excess energy to maintain T snow below 0°C. 3. Update T soil,n every 24 hr of the simulation (Equation 6). 4. Calculate d snow of the next time step using M snow and E s (Equation 1).

Determination of Model Parameters
Bulk ecosystem surface conductivity (k eco ), roughness lengths (z 0 and z T ), evaporation efficiency (β), and albedo (a) were evaluated as they represented canopy surface characteristics.These values were determined based on observation data, combining all data collected during the period from DOY 60 to SDD.The value of k eco was manually tuned for each year such that the modeled SDD agreed with the observed SDD.
Roughness lengths (z 0 and z T ) and the evaporation efficiency (β) were determined from observed values of friction velocity (u * ), H, λE, U, and T eco , and β.A single representative value for each of these parameters was determined for each site.First, the bulk transfer coefficients for surface stress (C M ) and heat (C H ) were estimated as the regression slopes using observed u * and H, along with their corresponding scalar gradients (e.g., Hayashi Water Resources Research 10. 1029/2023WR035984 et al., 2005)).Subsequently, roughness lengths and β were estimated using Equation 7. The observed value of T eco was estimated using observed upward longwave flux (R lu ) as Once C M and C H were determined from the regression, the values of z 0 and z T were determined by the following equations: The value of β was then calculated with the following equation: The average values of β across the years were determined for each site.
Albedo (a) was parameterized for each site based on the relationship between observed d snow and albedo (Appendix B).

Model Performance
We consider that the key to successfully modeling snowmelt is the accurate replication of d snow as well as each component of energy fluxes contributing to the determination of T snow (gray arrows in Figure 1).Modeled net shortwave flux (R s,net = (1 a)R sd ), εσT eco 4 , H, λE, and G 0 were compared with observation data as they determine T snow (Figure 1).Observed G 0 was calculated as the sum of observed ground heat flux and the heat storage estimated based on the uppermost T soil measurement.Although observations of λ f M snow and G s were not available, the accuracy of λ f M snow was inferred from d snow , and G s (=R s,net + εR ld εσT eco 4 -H λE) was considered proper if the other energy flux components were accurate.
The accuracy in terms of the magnitude of modeled variables was evaluated using the root mean square error (RMSE) for each site for both half-hourly and seasonal average values.The seasonal average values obtained from the model outputs were further assessed in comparison to the observation with respect to the accuracy in replicating interannual variations ("slope"), bias ("intercept"), and differences in performance between the sites ("site effect").This evaluation was conducted using the general linear model, using the glm function of R version 4.3.1 (R Core Team, 2023).

Sensitivity Analysis
To assess the effect sizes of meteorological drivers and of k eco on the interannual variations in SDD, the model was executed while maintaining one of the target variables at the 2011 values, as this year was considered representative of the average climatic conditions.This procedure is referred to as an experimental simulation, in contrast to the baseline simulation.
Subsequently, the difference in SDD between the experimental and baseline simulations (ΔSDD = Sensitivitydriven SDD-Baseline SDD) was calculated.A positive ΔSDD of a specific year indicates that a target variable of the year contributed to a delay of SDD.Finally, the absolute value of ΔSDD over the 10 years (mean |ΔSDD|) was computed to compare the effect sizes among different drivers.A greater value of the mean |ΔSDD| indicates a more significant impact of the target variable on snowmelt.
The evaluated meteorological drivers were categorized into two groups according to the model structure: halfhourly meteorological forcings (T air , q air , R sd , R ld , P snow , U, and atmospheric pressure) and initial conditions (d snow and T soil on DOY 60).The value of d snow on DOY 60 accumulated on the forest floor must be influenced not only by snowfall but also by canopy interceptions (Nakai et al., 1999).Nevertheless, we considered that the interannual variation in the initial value of d snow represented the winter snowfall, assuming the canopy structure did not considerably change over the decade.Likewise, we considered that the value of P snow represented spring snowfall.The sum of d snow on DOY 60 and total P snow represented the total snowfall on the forest floor for each season.The sensitivity analysis also tested the effect of the total snowfall on SDD.Furthermore, the effects of the differences in canopy surface parameters (k eco , roughness lengths, β, and albedo parameters) on the snowmelt timing were evaluated by replacing the parameters of US-Prr with the parameters tuned for US-Uaf.
Although different T air or q air could affect R ld (Ohmura, 2001;Rutter et al., 2023;Zhang et al., 1997), we treated these forcing variables as independent in this sensitivity analysis.The assumption would be reasonable if atmospheric radiation alone does not determine near-surface temperature conditions, as shown in the case of Sicart et al. (2008).

Snow Depth and Meteorological Conditions
Figure 2 shows the meteorological conditions in the spring of 2011, 2013, and 2016 as they represent the average, late, and early snowmelt.Table 3 lists the values of each observed variable for 10 seasons for both US-Prr and US-Uaf.Snowmelt occurred gradually, followed by a rapid disappearing phase for a few weeks (Figure 2a).The value of d snow at the beginning of March (DOY 60) and total P snow considerably varied among the years (Figures 2a,  Air temperature (T air ) was the second highest in 2016 after 2019 in both US-Prr and US-Uaf (Figure 2b, Table 3).The overall T air during the target period was lowest in 2013 in US-Prr and in 2017 in US-Uaf (Table 3).The value of q air increased as the snowmelt advanced (Figure 2c).Shortwave and longwave radiation (R sd and R ld ) conditions were similar between US-Prr and US-Uaf, but R ld varied more than R sd between years (Figures 2d and 2e).The average R ld was greatest in 2016 in US-Prr.Albedo (a) reflected the trajectory of d snow (Figure 2f).SDDs were 5 days on average earlier in US-Uaf than in US-Prr, except that the difference was not discernible in 2016 (Table 3).

Canopy Surface Properties
The determined k eco values ranged from 0.37 W m 1 K 1 in 2016 to 3.0 W m 1 K 1 in 2012 at US-Prr and from 0.5 W m 1 K 1 in 2016 to 8.2 W m 1 K 1 in 2018 at US-Uaf (Table 4).The values of k eco tended to be smaller in the years with less snow (initial d snow and P snow in Table 3) (e.g., 2016 and 2019) compared with the years with relatively greater amounts of snow (e.g., 2018 and 2020).The calculated z 0 and z T based on Equation 10were greater for US-Uaf than the values determined for US-Prr (Table 2) as the differences between the sites were obvious in C M and C H represented as the slope of the regression lines in Figures 3a and 3b.The evaporation efficiency (β) was greater in US-Prr than in US-Uaf (Table 2).

Model Performance
The model was run using the determined model parameters, and the results were compared with observed data (Figure 4).The value of k eco was tuned to agree modeled and observed SDD (Figure 4a).The results for 2013 are shown as an example in Figure 4  , and G 0 (Figures 5a-5c, f, Table 5).The model accurately represented the relative magnitudes of H and λE, with their respective RMSEs of 11 W m 2 and 4 W m 2 , which were not notably high when compared to the RMSEs of other energy flux components (Figure 5).While the interannual variations in H and λE did not precisely match the observations (as seen in the slope effect in Table 5), both the modeled and observed H consistently displayed their highest values in 2013 and low values in 2019.Note that H data in 2019 from US-Uaf were not utilized due to large gaps in the data.

Sensitivity Analysis
We found the effects of snowfall (i.e., initial d snow and/or P snow ), T air , and R ld were key meteorological drivers that determine the interannual variations in SDD at US-Prr (Figure 6, Table 6).The initial d snow at the beginning of the model simulation (DOY 60) was identified to be the most important driver (mean |ΔSDD| = 9 days).The total snowfall effect, combining effects of P snow and initial d snow , accounted for 10 days of the interannual variation in SDD.Other than the factors related to snowfall, T air was also an important meteorological driver (mean | ΔSDD| = 7 days) followed by the effects of R ld (mean |ΔSDD| = 5 days).Humidity, solar radiation, wind speed, and atmospheric pressure effects were minimal even if combined (mean |ΔSDD| = 2 days) ("Other met.forcings" in Figure 6 and Table 6), suggesting that the interannual variations in these variables had minimal influence on  Important meteorological drivers were different by extremely early and late snowmelt years.The role of T air relative to other meteorological drivers was apparent in 2013 when the low T air resulted in the delay of the simulated SDD by 11 days (Figure 6, Table 6).While the year 2013 had a relatively large amount of spring snowfall (i.e., P snow ) (Table 3), the total snowfall effect (i.e., combined P snow and initial d snow effects)  3).(ΔSDD = 5 days) was smaller than the T air effect.On the other hand, the low amount of the accumulated winter snow (i.e., initial d snow ) was important for explaining early snowmelt in 2016 and 2019 accounting for 3 weeks of the advancement of SDD (Figure 6, Table 6).Although the effect of R ld on the snowmelt was minimal in 2013, the high R ld advanced SDD by 11 and 13 days in 2016 and 2019, respectively (Table 6).

Water Resources Research
The model also showed that the variations in k eco and the parameters for albedo were as crucial as meteorological factors in determining snowmelt timing (Figure 7, Table 6).Simulations using the fixed k eco of 2011 resulted in an average error of 8 days in SDD (Figure 7a, Table 6).Simulations altering only k eco from US-Prr to US-Uaf led to an estimated difference in SDD of 7 days (Figure 7b, Table 6).Using the albedo parameters of US-Uaf advanced the snowmelt by a week (Table 6).The effects of roughness lengths and evaporation efficiency were minimal (ΔSDD was 3 days or less).IKAWA ET AL.
To further evaluate T air effects on snowmelt, the differences in modeled H and H of the simulations with fixed T air in 2011 were compared (Figure 8).The simulation suggested that a lower level of T air tended to increase H cooling the snow surface, which was particularly apparent in 2013.On the other hand, higher levels of T air tended to decrease H, enhancing snowmelt in the years 2016 and 2019.

Snowfall or Other Spring Meteorology?
To our knowledge, the information is limited with respect to the importance of snow depth or snowfall compared with other meteorological drivers for the snowmelt timings.This could be because many of the previous studies had focused on comparing the magnitudes of different energy fluxes utilized for snowmelt (e.g., Hayashi et al., 2005;Langer et al., 2011;Nishimura et al., 2018;Pomeroy et al., 2009;Suzuki et al., 2006;Zhang et al., 1997).The decadal observation data that we presented in this study include the information on the naturally occurring interannual variations in snowfall and heat fluxes.Therefore, applying the snowmelt model to the decadal data sets allowed us to investigate not only the distributions of the energy fluxes illustrated in Figure 1 but also the relative importance of snowfall and other meteorological drivers.
Our results indicated that the winter snowfall (i.e., initial d snow ) was the most important driver explaining 9 days of the interannual variations in SDD.Together with the spring snowfall (i.e., P snow ), the total amount of snowfall explained 10 days of the interannual variations in SDD.According to Ballinger et al. (2023), there has been a notable increase in air temperature and reduced spring snowfall observed from 1957 to 2021.These changes would have contributed to the advancement of spring snowmelt.They also reported an increasing trend of winter snow water equivalent, which offsets the advancement to a certain extent.As the flux observations from several  sites are expected to span multiple decades in the near future, the presented modeling framework is anticipated to facilitate an attribution analysis covering this extended time period.
Although our model-based approach identified that the snowfall had the greatest effect on predicting the SDD overall, spring meteorological conditions such as T air and R ld were also important for the snowmelt processes.For example, the low T air characterized the late snowmelt in 2013.Even though the year 2013 experienced a high amount of P snow , the total snowfall effect (ΔSDD = 5 days) was less than the low T air effect (ΔSDD = 11 days) (Table 6).

Air Temperature or Longwave Radiation?
Air temperature and atmospheric radiation values are often highly correlated, and it is difficult to separate their individual impacts on snowmelt (Ohmura, 2001).Many studies in high-latitude ecosystems have documented the importance of longwave radiation in the snowmelt process (Suzuki et al., 2006;Suzuki & Ohta, 2003;Zhang et al., 1997).Zhang et al. (1997) reported that R ld was a key meteorological driver for the snowmelt in the Arctic and Subarctic based on the results of an atmospheric radiative transfer model.Water Resources Research 10.1029/2023WR035984 Our results indicated that the H effect via the change in T air had a greater effect than the R ld effect on the interannual variations in the snowmelt (Figure 6, Table 6).While the importance of sensible and latent heat fluxes in the snowmelt process has been reported (Hayashi et al., 2005;Langer et al., 2011;Nishimura et al., 2018;Sicart et al., 2008), little is known about the extent to which these turbulent fluxes affect the interannual variation in snowmelt timings.In our study, the cool spring in 2013 and the warm spring in 2019 were found to characterize high or low values of H, respectively (Table 3, Figures 5, 6, and 8).This finding resembles Sicart et al. (2008), who reported that high variability in temperature caused variations in the glacier melt through sensible heat flux.High T a conditions are often associated with high R sd or R ld conditions, which enhance H, making it challenging to isolate T a effects solely through observations, as exemplified in the warm spring of 2016.The model-based approach exemplified here is effective in quantitatively isolating the T a effect on H, as depicted in Figure 8.We also consider that the importance of H emerged in this study because of the high z 0 compared with a flat snow surface (e.g., 1.4*10 4 in Kondo, 1994) and the much greater sensible heat transfer compared to the latent heat transfer as represented in the value of β.
We considered the primary effect of different T air was manifested through different H.However, it is important to note that different T air affects snowmelt through other factors, resulting in diverse effects.For example, despite the relatively low T air in 2017, its effect on SDD was positive (ΔSDD = +4 days) compared with 2011, indicating the advancement of snowmelt.Our sensitivity analysis suggested that lowering T air tended to decrease the loss of the longwave flux from the surface (εσT eco 4 ) as T eco decreased, thereby enhancing snowmelt (Text S5 in Supporting Information S1).Using different T air values also influenced snowmelt through changes in λE, but this effect was relatively smaller compared to the effects of H or longwave radiations.
It should also be noted that an important environmental process for snowmelt is also different according to a research question to be asked.Suzuki et al. (2006) identified the importance of net radiation in the mountainous region of eastern Siberia when comparing snow ablation between a larch forest and an open field.Our results also suggested that atmospheric radiation had greater effects than the turbulent heat flux when solely focusing on the magnitude of the energy fluxes, rather than on their interannual variations.Several studies suggested the importance of longwave radiation within forest canopies (Hashimoto et al., 1998;Pomeroy et al., 2009;Rutter et al., 2023;Suzuki & Ohta, 2003;Suzuki et al., 1999).While our study addresses solely on integrated ecosystem fluxes, the distribution of thermal radiation and turbulent heat fluxes may differ on the snow surface.This difference could lead to different interpretations regarding the distributions among radiation and turbulent fluxes when assessed in proximity to the snow surface.

Extremely Cold and Warm Springs
While the snowfall, T air , and R ld were identified as key meteorological drivers related to the interannual variations in the snowmelt timings, their relative importance varies from year to year.This section focused on the cold spring of 2013 and the warm spring of 2016.
The sensitivity analysis suggested that the late SDD in 2013 was attributed to low T air , which increased H (Figure 8), cooling the snow surface.Relatively high spring snowfall (see P snow in Table 3) in 2013 further contributed to the delay of snowmelt (Figure 6 and Table 6), but the total snowfall effect (ΔSDD = 5 days) was smaller than the low T air effect (ΔSDD = 11 days) according to Table 6.This result underscores that the relative importance of each meteorological driver on the snowmelt timings depends on the meteorological conditions in spring.
At the continental scale, the extremely warm winter in 2016 was attributed to the sparse snowpack that decreased the albedo (Walsh et al., 2017) and the greater downwelling longwave radiation (Cullather et al., 2016).According to the near-ground evaluation in this study, the low amount of the winter snowfall (i.e., initial d snow ) was attributed to 3 weeks (ΔSDD = 22 days) of the advancement of SDD.Walsh et al. (2017) reported the statewide average snowfall in February 2016 was 20% of the normal year, which was consistent with our observation of the low initial of d snow in 2016.While the effect T air was still important for the early snowmelt in 2016 (ΔSDD = 11 days), the role of R ld was as important as T air (ΔSDD = 11 days).The emerging importance of R ld was also confirmed in another warm year, 2019.

Usability and Limitation of the Bulk-Surface Approach
Introducing the concept of a bulk surface has enabled us to understand important processes related to snowmelt within a simple energy balance framework, which considers only seven energy fluxes to be resolved: R s,net , εσT eco 4 , H, λE, G 0 , G s , and λ f M snow .All in-canopy processes were aggregated into a few parameters, namely albedo parameters, roughness lengths, β, and k eco , which can be parameterized based on ordinary eddy covariance site data.This approach can bypass the requirement for detailed canopy structure information, which is often challenging to acquire in various forest sites, while still maintaining theoretical correctness from the perspective of energy balance.Despite the simplicity, the model successfully replicated d snow as well as the magnitudes and variations of each energy flux component.We further confirmed that our model essentially produced results similar to those of the Flexible Snow Model (FSM2: https://github.com/RichardEssery/FSM2)(Mazzotti et al., 2020), which incorporates explicit canopy parameterizations.Specifically, both low T air and high R ld clearly affected the temporal patterns of d swe and SDD (Text S6 in Supporting Information S1).
It is important to note that the primary purpose of the bulk-surface model is to interpret the underlying processes of observed data within a theoretical framework, rather than serving as a general-purpose model.Hence, we did not model intercepted snow, which modifies energy balance processes (Lundquist et al., 2021;Suzuki & Nakai, 2008).Instead, we utilized already observed snow depth data (i.e., d snow ) to estimate the snowfall rate on the forest floor (P snow ), focusing only on the period when the snow accumulated on the canopy was minimal in spring.The inclusion of a canopy interception sub-model (e.g., Pomeroy et al., 1998) is necessary if the model is intended for the extended period from snow accumulation or in locations where there is no available information about snowfall on the forest floor.
Not considering the canopy and forest floor led to some issues requiring close attention upon implementing the bulk-surface approach.The primary issue is an inevitable drawback: the bulk surface approach does not differentiate between different temperature conditions among the canopy, forest floor, and snow.This likely resulted in a disagreement of modeled and observed H when the temperature between the canopy and snow surface would deviate near SDD.Considering at least two distinct canopy compartments, as Mazzotti et al. (2020) demonstrated, would enhance the model's accuracy, though it would necessitate additional canopy information.Flux-partitioning methods (e.g., Constantin et al., 1999) together with the lidar application (Mazzotti et al., 2020;Pennypacker & Baldocchi, 2016) might give us a hint to achieve parameterizations based on eddy covariance data for a process-based model with explicit canopy consideration, such as the Flexible Snow Model (Magnusson et al., 2019;Mazzotti et al., 2020).
The second issue is that we had to introduce an empirical parameter, k eco , for which we lack a clear physical meaning, yet it has a significant impact on the estimate of SDD.The parameter of k eco in our model would partly explain the ablation effect, as its value tends to be greater in the years with greater amount of snowfall.The socalled decoupling between above and within the canopy often occurs in boreal forests (Rutter et al., 2023;Starkenburg et al., 2013).It is reasonable to consider that such an effect is also encompassed within the value of k eco .However, the lack of a clear physical interpretation of k eco limits the utility of our model for evaluating incanopy processes and for generalizing the model to make it readily applicable to various other sites.Despite the uncertainty associated with k eco , the sensitivity analysis using a fixed k eco value in 2011 across all years was still useful for qualitatively estimating the change in SDD (ΔSDD) in response to environmental stimuli at our sites.
The third issue is that we utilized a simplified snow scheme that utilizes a constant albedo throughout a day.Transmission and reflections on the snow surface are influenced by solar zenith angle (Aoki et al., 2000).Our model also treats the snow as a single layer, even though the morphology of snow varies with depth (Niwano et al., 2012;Saito et al., 2012).Considering the detailed albedo and morphology of snow would be necessary if factors affecting the fate of shortwave radiation at a sub-daily timescale are relevant to the main research questions.

Conclusions
A simple snowmelt model was developed to interpret the interannual variations in the spring snowmelt timings in the open black spruce forests in interior Alaska.The model-based approach suggested that snowfall, T air , and R ld were key meteorological drivers of the interannual variation in spring snowmelt in the black spruce forests in Alaska.The effects of k eco and the parameters for determining albedo were also important determinants of spring snowmelt timings.The magnitude of H was relatively small in the energy budgets, but the increase in H under low T air in 2013 considerably delayed SDD.The bulk-surface approach presented here does not resolve the processes occurring in proximity to the snow surface.Nevertheless, we have demonstrated that this approach still offers valuable insights into important processes related to observed snowmelt data within a simple energy balance framework.
where a max and a min are the maximum and minimum albedos before and after the snowmelt period, respectively, and γ 1 , γ 2 , γ 3 are fitting constants.The value of a min was fixed at 0.3 corresponding to the criteria we determined SDD.The values of γ 1 , γ 2 , γ 3 , and a max were determined using the nonlinear regression analysis (nls function of R).The value of a was generally greater in US-Prr than in US-Uaf at a given snow depth during the rapid snowmelt period (Figure B1).

Figure 1 .
Figure 1.Schematic description of the snowmelt model (a) and the calculation steps (b).

Figure 2 .
Figure 2. Observed half-hourly snow depth (d snow ) (a), and daily values of observed air temperature T air (b), specific humidity q air (c), solar radiation R sd (d), atmospheric radiation R ld (e), and albedo a (f) for March-May (DOY 60-151) in 2011, 2013, and 2016 in US-Prr.Vertical lines indicate snow disappearance date for each year.

Figure 3 .
Figure 3.The relationships of observed friction velocity u * (a) and sensible heat flux H (b) and their driving variables under the near-neutral condition (|z obs /L| < 0.1, where L is the Obukhov length) for DOY 60-SDD.Wind speed and air temperature were corrected for the height of 10 m (U 10 and T air,10 , respectively) to allow the direct comparisons between the two sites.

Figure 4 .
Figure 4. Observed and modeled values of (a) snow depth (d snow ), (b) net shortwave radiation (R s,net ), (c) upward longwave radiation (εσT eco 4 ), (d) sensible heat flux (H), (e) latent heat flux (λE), and (f) the heat flux from snow to soil (G 0 ) in 2013 in US-Prr.The year 2013 data were selected as examples due to the extensive duration of the snow cover (Table3).

Figure 6 .
Figure 6.Comparisons between the sensitivity-driven snow disappearance date (SDD), where individual variables were set at the value of 2011, and the baseline SDD driven from the baseline simulation at US-Prr for (a) meteorological forcings, (b) model initial conditions, and (c) "Total snowfall," representing the total snowfall in the season, combining effects of "Spring snowfall" and "Initial snow depth."The departure from the 1:1 line indicates the effect size of each variable."Other met.forcings" represents the combined effects of specific humidity, solar radiation, wind speed, and atmospheric pressure.

Figure 7 .
Figure 7. Comparisons of the sensitivity-driven snow disappearance date (SDD) and the baseline SDD driven from the baseline simulation for US-Prr (a) when k eco was fixed at the value of 2011 to evaluate the effect of the interannual variation in k eco on SDD and (b) when k eco of US-Uaf was used instead to evaluate the effect of the site difference in k eco on SDD.

Figure 8 .
Figure 8.The differences between modeled sensible heat fluxes (H) and the H when the input T air was fixed at the values of 2011 [H(T air,2011 )] compared with the mean T air at US-Prr.The data in 2011 is marked with the white, indicating it is the reference.The output values were averaged over 41 days (DOY 60-100).

Figure B1 .
Figure B1.The relationships between the daily average albedo and the snow depth (d snow ) at US-Prr in 2011-2020 and US-Uaf in 2017-2020.

Table 1
Information of Instruments

Table 2
Variables and Constants Used for the Snowmelt Model IKAWA ET AL.

Table 3 )
. SDD was the earliest in 2016 (DOY 98) coinciding with the lowest snowfall (d snow on DOY 60 + Total P snow = 0.35 m) for US-Prr.The earliest SDD was in 2019 followed by 2016 in US-Uaf.SDD was the latest in 2013 for both US-Prr (DOY 133) and US-Uaf (DOY 131).

Table 3
Summary of Observation Data Year T air q air R sd R ld H λE d snow on DOY 60 Total P snow d snow on DOY 60 +Total P snow SDD Note.Average values are listed for T air (°C), q air (g kg 1 ), R sd (W m 2 ), R ld (W m 2 ), H (W m 2 ), and λE (W m 2 ) for DOY 60-100.The period from DOY 60 to DOY 100 effectively characterizes surface flux and meteorological conditions during each snowmelt season.The value of d snow on DOY 60 represents snowfall in winter, the total P snow (m) represents total snowfall in March-May (DOY 60-151), and their sum (d snow on DOY 60 + Total P snow ) represents the total snowfall in the season.IKAWA ET AL.
because the year 2013 had the longest simulation period from DOY 60 until SDD.Overall, we considered the model was able to reasonably capture the temporal trajectories of observed d snow , R s,net , εσT eco 4 , H, λE, and G 0 and their temporal variations given the ranges of the observed values.The root mean square error, RMSE, for d snow averaged for all years combining both sites was 0.06 m.Likewise, the RMSE of R s,net , εσT eco 4 , H, λE, and G 0 were 15 W m 2 , 13 W m 2 , 33 W m 2 , 9 W m 2 , and 17 W m 2 .RMSEs of individual year and site are listed in Text S3 in Supporting Information S1.There was a noticeable disagreement between the model and observations for εσT eco 4 and H near SDD.The overestimation of the diurnal amplitude of G 0 was mitigated if a lower k soil value was assigned to the topmost soil layer (data not shown).Still, we chose not to do so to maintain the simplicity of the model.The model captured the magnitude and interannual variations of the seasonal average value of d snow , R s,net , εσT eco

Table 4
Bulk Ecosystem Surface Conductivity (k eco )

Table 5
The Root Mean Square Error (RMSE) and the General Linear Model Test for Seasonal Average Values of d snow , R s,net , εσT eco 4 , H, λE, and G 0 for the Combined Data From US-Prr and US-Uaf

Table 6
Departures of the Sensitivity-Driven Snow Disappearance Date (SDD) From the Baseline SDD (ΔSDD) and Mean |ΔSDD| for US-Prr.