Global Simulation of Snow Algal Blooming by Coupling a Land Surface and Newly Developed Snow Algae Models

Snow algae are found from spring to summer on snowfields and glaciers throughout the world. Their blooming darkens snow surfaces, reducing snow surface albedo and accelerating melting. Uncertainties remain, however, regarding the blooming season and global distribution of these algae. To reproduce snow algal bloom temporal and geographical variability, we improved an existing snow algae model using a land surface model calibrated with a reanalysis dataset of the global atmosphere. Snowfall and daylight length data for selected model locations were also incorporated. To evaluate its performance, we used in situ observational data from 15 polar to alpine area sites. The improvements made in this study allowed the reconstruction of detailed snow algal blooming reports from various locations worldwide, and the results suggested that the major factors affecting the appearance of snow algal blooming were the snow melting period duration and algal growth interruption by new snow cover. We then incorporated the updated snow algae model into a land surface model and performed a global simulation. In this case, our simulation suggested that red snow could appear on snowfields during the melting season but only in the absence of frequent new snowfalls, and if the snow cover persists long enough to allow prolonged algal growth. This simulation has the potential to be used for global prediction of future red snow phenomena, which are likely to synchronize with global climate change.

A number of land surface models, which can be driven using reanalysis data of the global atmosphere, have been proposed to simulate temporal and spatial changes in snow properties on a global scale. For example, CLM (Lawrence et al., 2019), ORCHIDEE (Krinner et al., 2005), ISBA (Decharme et al., 2016), and Minimal Advanced Treatments of Surface Interaction and Runoff land surface model (MATSIRO; Nitta et al., 2014Nitta et al., , 2017Takata et al., 2003) have been established as land surface models incorporated into climate models to represent physical land processes. These models can calculate temporal and spatial changes in the snow water equivalent, snow  (Cepák et al., 2016;Fujii et al., 2010;Ganey et al., 2017;Hisakawa et al., 2015;Huovinen et al., 2018;Lutz et al., 2014Lutz et al., , 2015Moestrup et al., 2018;Novis, 2002;Onuma et al., 2018;Painter et al., 2001;Remias et al., 2016;Spijkerman et al., 2012;Stibal et al., 2007;Takeuchi & Kohshima, 2004;Tanaka, 2016;Tanaka et al., 2016;Thomas & Duval, 1995;Vimercati et al., 2019). Solid, open, and open-dotted marks indicate seasonal data validation sites, 1-day data validation sites, and no applicable sites in this study, respectively. Only validation sites where red snow algal abundance was quantified by the cell count method were used. Site names correspond to those in Table 1. 10.1029/2021JG006339 3 of 22 temperature, water runoff, evaporation, and sublimation globally. In addition, a simulation with a land surface model can be conducted independently, using atmospheric conditions near the land surface as the input data. Such a simulation is termed an offline land simulation. Reanalysis data are obtained by simulations of the atmosphere to be consistent with meteorological observation, which creates gridded historical data sets at global or regional scales. The data are generally used as atmospheric conditions for offline land simulations, which allows temporal and spatial changes in land physical properties to be reproduced without model bias that may derive from atmospheric conditions.
In this study, we first improved the snow algae model using in situ snow algal abundance data, as reported from 15 locations worldwide, and the physical and meteorological snow conditions for these locations, as obtained from a land surface model. We incorporated the effect of snowfall and daylight length into the snow algal model established by Onuma et al. (2018) and then conducted offline land simulations at the study sites using various atmospheric conditions and biological parameters. Finally, we performed a global simulation of the land surface model, including the snow algae model, using reanalysis datasets of the global atmosphere to investigate seasonal and geographical variations in snow algal blooms worldwide.

Snow Algae Model
We used a snow algae model to calculate temporal changes in the abundance of snow algae across various snowfields. Temporal changes in the abundance of S. nivaloides on surface snow can be expressed using a differential logistic growth equation. The population density of the microbes was calculated as shown in Equations 1 and 2 (Onuma et al., 2018): where X represents the population density of microbes at time t, and μ indicates the hourly microbe growth rate. K denotes the snow surface algae carrying capacity. In this model, X increases when the snow surface temperature, T sn , is above 0°C (as algal growth only occurs on the melting snow surface). When algae first appear on the snow surface, X represents the initial population density X 0 . Although algal cells observed in the red snow surface are often in the cyst stage (e.g., Onuma et al., 2018)-when their populations do not actively increase-the model assumes algal growth on the snow surface, including the condensation of algal cells grown at the subsurface, with snow melt.
In this study, we improved this model, which may include growth and/or condensation of the algal cells, to broadly reconstruct in situ observations of algal cell abundance reported for snow surfaces worldwide. To calculate temporal changes in the abundance of red snow algae on surface snow at various locations, we added the effects of snowfall and daylight length into the original snow algae model proposed by Onuma et al. (2018). Onuma et al. (2016) reported that the abundance of a snow alga Chloromonas (C.) nivalis on a snowpack in Japan decreased when there were occasional snowfalls in spring; however, the snow algae model was not able to simulate such a decrease because it assumed a monotonic algal abundance increase. Because the snow algae model could calculate temporal changes in algal cell concentration in a surface snow layer to a depth of 2 cm, the accumulation of new snow above the algal layer should result in a decline in algal cell concentration. In this study, we updated the snow algae model to quantify the effect of snowfall on the algal abundance in the top 2 cm layer of surface snow.
Daylight length is another metric with the potential of affecting algal growth because snow algae grow photosynthetically. A previous study reported that snow algal blooms first appeared under light conditions with a penetration of 0.1% of the surface radiation (Curl et al., 1972), suggesting that snow algae on snow surfaces 10.1029/2021JG006339 4 of 22 can grow during daylight but not during the night. Therefore, we assumed that snow algae grow in sunlight and incorporated a day-night cycle effect on snow algal growth into the snow algae model, as shown in Equation 3: where X can increase when the snow surface temperature, T sn , and downward shortwave radiation, SW, are above 0°C and 0 W m −2 , respectively, at t.
We added two further equations to the calculation of snow algal growth, to quantify the effect of snowfall on snow algal abundance, as shown in Equations 4 and 5: where f dp indicates decreasing term as a function of the snowfall rate PR sn (mm s −1 ). In this model, we assumed that a new snow layer, which not include algal cells, is added to the top of the snowpack when snowfall occurs. Because the model calculates the algal cell concentration at the top 2 cm layer of surface snow, X in the top 2 cm should be reduced with an increase in the depth of the new snow layer. Equations 4 and 5 calculate X reduced by new snow cover in this study. After then, an increase in the algal abundance is calculated using Equations 1-3 based on the reduced X. If the new snow cover, which is calculated using PR sn and snow density D sn (kg m −3 ), exceeds 20 mm in 1 hr, X is reset to X 0 . In this study, D sn was assumed to be 300 kg m −3 , which is generally used for this parameter in the land surface model MATSIRO used in this study (Nitta et al., 2014;Takata et al., 2003). The land surface model is described in Section 2.3.
An overview of the original and two updated snow algae models can be seen in Figure 2. In this study, we evaluate the performance using the three snow algae models: Ⅰ original snow algae model shown in Equations 1 and 2 (X M ), Ⅱ the former updated model shown in Equations 1-3 (X MD ) and Ⅲ the later updated model shown in Equations 1-5 (X MDF ). The former updated model is a snow algae model to simulate algal cell concentration including the effects of snow melting and daylight length. The latter updated model is the former plus the effect of snowfall. In addition, we defined algal growth period GP to quantify the duration required for the specific algal cell concentration. The GP is referred to as GP M , GP MD and GP MDF for X M , X MD and X MDF , respectively. For example, GP MD represents the cumulative hours when T sn is above 0°C under daylight conditions. In the case of GP MDF , the cumulative time is obtained from X MDF by back calculating the solution of Equation 1. In the X MDF model, the initial cell concentration, X 0 , becomes established at the first snow appearance, and then the cell concentration increases with the duration of snow melting under daylight conditions before decreasing when new snowfall covers the surface ( Figure 2). As the model could not consider the migration of motile algal cells in a snowpack, algal cell concentrations reduced by snowfall were not conserved in the model snow layers and did not increase again until there was further snow melting under daylight conditions. The modeled algal cell concentrations reverted to zero once the snowpack disappeared.

Using In Situ Data for Snow Algal Abundance to Validate the Snow Algae Model
In this study, we used temporal changes in algal cell abundance observed in situ at 15 locations as model validation data. An overview of the validation sites located in polar and alpine regions-six sites representing sites used in previous studies, and nine representing those observed specifically for this work-can be seen in Figure 1c and Table 1.
Surface snow algal abundance data reported from Svalbard (Arctic Ocean; SV-W and SV-S; Lutz et al., 2014;Stibal et al., 2007), SE Greenland (GL-A, Lutz et al., 2014), Europe (EU-T, Remias et al., 2016), North America (NA-S, Painter et al., 2001), and New Zealand (NZ-A, Novis et al., 2002) were also used. Because most algal cell abundance data had been reported as algal cell concentration per unit of melt water volume (cells L −1 ), we converted these data to algal cell concentrations per unit of area, assuming that the snow density was assumed to be 500 kg m −3 and that the collected samples had been 2 cm deep. The assumed value is closed to snow density when red snow occurs reported by previous studies (Müller et al., 2001;Onuma et al., 2018).
We also used unpublished algal data obtained from snowfields in Japan (Mt. Tateyema) and the Antarctic Peninsula (Livingston Island), in 2012 and 2015 respectively. Mt. Tateyama (N 36.3°, E 137.4°) is an alpine snowfield above the tree line (1,850-3,000 m ASL), located in W Japan, and red snow algal blooming can be observed there annually as the snow pack thaws (May-July; Segawa et al., 2005). We selected a snowy plateau called Raichozawa, located at 2,300 m, as the study site (JP-T) because the snowpack here remained until early August.  Island from January-February . Our study site (AP-L) was located on the Livingstone Island coastal snowfield.
Temporal changes in the abundance of red snow algae were quantified using surface snow samples collected through days 117-217, in 2012, at JP-T, and during 6 sampling days during the period extending through days 19-30, in 2015, at AP-L. Samples were collected on each observation date from 3 to 10 randomly selected surface locations (0-2 cm depths), using a stainless-steel scoop. The sampling areas occupied ∼100 cm 2 and were recorded for each collection. Samples were melted on-site and preserved in 3% formalin in 30 mL clean polyethylene bottles before being transported to Chiba University, Japan, for analysis. Algal abundances were obtained from the water samples by counting cells and were represented as the cell number per unit surface area of snowpack (cells m −2 ). This methodology has been described in more detail in previous studies, which included those in which field observations were carried out at sites GL-Q, AK-E, SB-S, AK-K, AT-A, TN-U, and PG-T, as shown in Table 1.

MATSIRO Land Surface Model
We used the MATSIRO (Nitta et al., 2014(Nitta et al., , 2017Takata et al., 2003) to evaluate temporal changes in snow algal abundance on snowfields worldwide. This model was developed to simulate land-based physical processes in a general circulation model, and six versions-the MATSIRO6, involving up to three snow layers, six soil layers (14 m in total), and a single canopy layer-were used for the Model for Interdisciplinary Research on Climate (MIROC6; Tatebe et al., 2019).
MATSIRO6 can simulate temporal and spatial changes in the snow water equivalent, snow cover fraction, and snow temperature. The snow water equivalent was simulated based on the water balance, and in MATSIRO6 it was derived from the snowfall rate, snow sublimation, snowmelt, and refreezing of rainfall and snowmelt. The snow cover fraction was simulated using a lognormal distribution function for the subgrid snow water equivalent distribution, whereas the temperature of each snow layer was simulated using a thermal conductivity equation. The detailed methodology for calculating these physical properties of snow may be found in Nitta et al. (2014). Temporal changes in these snow physical properties have been validated using observations from various snowfields and derived from data for other land surface and snow physical models reported by the model intercomparison project (ESM-SnowMIP; Krinner et al., 2018). In this study, we used snow surface temperatures calculated using MATSIRO6 as input data for snow algal simulations at the study sites.

Atmospheric Conditions for Land Surface Modeling
The atmospheric conditions used in MATSIRO6 simulations were derived from reanalysis data ( Table 2). Various datasets for reanalysis-which were derived from reanalysis data of global atmosphere near the land surface, and bias-corrected using global meteorological observations-have been established for land surface modeling. In this study, we used the atmospheric conditions derived from the reanalysis dataset for each study site, because time-series meteorological observations were not always available.
The WFDEI reanalysis dataset of the global atmosphere (Weedon et al., 2014) was used as the atmospheric conditions for the land surface modeling in this study. This reanalysis dataset includes three-hourly information on surface air temperature, surface air pressure, downward radiation (shortwave and longwave), humidity, wind speed, and precipitation rate. While this reanalysis dataset is appropriate for land surface modeling at global or regional scales, high levels of uncertainty are present when it is used for meteorological conditions at specific elevations due to its rough horizontal resolution (0.5° × 0.5° globally).
For this reason, we applied elevation corrections to simulate temporal changes in snow algal abundance at specific sites. The surface air temperature (T atm ) at each site was corrected using elevation information and the original air temperature and by applying a temperature lapse rate, which was assumed to be −6.5 × 10 −3 K m −1 ( Figure S1 in Supporting Information S1). The surface air pressure (PS atm ) of each snowfield was corrected using elevation information and surface air temperature (before and after the correction), as recommended by the World Meteorological Organization (WMO-No. 8, in CIMO Guide, Part I, Chapter 3). The specific humidity (Q atm ) at each site was corrected from the original data using the ratio of surface air pressure before and after correction, whereas the snowfall rate (PR sn ) was corrected from total precipitation data (PR), using the ratio of rain to snow, which itself was estimated using the surface air temperature, air pressure, and specific humidity given in MATSIRO6. The input and output data for the offline land simulation used to calibrate MATSIRO6 in this study can be seen in Table 2

Experimental Design of Snow Algal Growth Simulation at 15 Different Locations
To evaluate the algal model version as improved in this study, algal growth simulations (Ag-exp) were conducted under three different conditions at 15 sites, as shown in Figure 1c and Table 1. The three conditions covered snow algal growth with (a) the effect of snow melting only (X M ), (b) effects of snow melting and daylight length only (X MD ), and (c) effects of snow melting, daylight length, and snowfall (X MDF ; Table 3). The atmospheric conditions for the land offline simulations were supplied by the WFDEI reanalysis dataset, which had been corrected for study site elevations in advance.
The initial snow depth (in water equivalent) in MATSIRO6 was assumed at each study site to remain stable until the maximum algal cell concentration date. For example, the initial value was assumed to be 10,00 kg m −2 for Arctic sites and 3,000 kg m −2 for Japanese sites, which was consistent with snow depth observations (equal to 9 m in winter at JP-T, Osada et al., 2004). The snow surface temperature (T sn ) calculated using MATSIRO6 and the derived snowfall rate (PR sn ) and solar radiation (SW) from reanalysis data were used as input data for the snow algae model in this study (Figures S2 and S3 in Supporting Information S1).
Biological parameter data-such as initial cell concentration (X 0 ), growth rate (μ), and carrying capacity (K)were not generally available for red snow algae in snowfields worldwide. To overcome this, we applied field observation values from the snowfield of a Greenlandic glacier, as reported by Onuma et al. (2018), for Ag-expand so X 0 and μ were assumed to be 6.33 cells m −2 and 0.018 hr −1 , respectively. Maximum algal cell concentrations observed from the sites (Table 1) were used for K because the carrying capacity of red snow algae might vary for each snowfield, as suggested by Onuma et al. (2018).
Simulations were conducted from January 1 to December 31 each year at the study sites, except for the southern hemisphere sites NZ-A and PG-T, where 1998 and 1999 data for July 1 to December 31 were used.
We conducted more sensitivity testing on the updated snow algae model using two biological parameters-initial cell concentration (X 0 -exp), and growth rate (μ-exp)-at each observation site (Table 3). Onuma et al. (2018) reported that the S. nivaloides initial cell concentration and growth rate were 6.33 and 694 cells m −2 and 0.39 and 0.42 days −1 (0.016 and 0.018 hr −1 ), respectively, at 2 sites on a Greenlandic glacier. Field observations suggested that S. nivaloides algal spores (cysts) were wind induced onto the snow surface during the early melting season, and they assumed that the initial cell concentration consisted of these wind-supplied algal spores before on-site algal growth initiation. Based on a previous study, the range for the initial cell concentrations (X 0-exp) was assumed for sensitivity testing to be between 1-1,000 cell m −2 . Because snow algae growth rates vary between snowfields worldwide, it had been suggested previously that the sensitivity of the snow algae model to growth rates should be investigated. Previous studies have reported that C. nivalis growth rates, as obtained from field observations and cultivation, were 0.22 (Onuma et al., 2016) and 0.60 days −1 (Leya et al., 2009). Based on these data, we assumed that the μ-exp growth rate ranged between 0.01-0.025 h −1 . The maximum algal cell concentration observed from the study sites was assumed to be K, for both X 0 -exp and μ-exp, the same as for Ag-exp.

of Sensitivity Tests Applied in This Study
and CRUJRA reanalysis datasets, was conducted in this study, to review the sensitivity of algal growth to atmospheric conditions (Table 3). GSWP3-FD (Hurk et al., 2016;Kim, 2017) and CRUJRA (Harris, 2019) reanalysis datasets provided the different atmospheric conditions, and an overview of the three datasets, including WFDEI, can be seen in Table S1 in Supporting Information S1. The surface air temperature, surface air pressure, and specific humidity data from GSWP3-FD and CRUJRA were corrected using each site's elevation, as was WF-DEI. Besides, we additionally conducted Lr-exp, which is an ensemble simulation using temperature lapse rate, ranging from −8.0 × 10 −3 to −5.0 × 10 −3 K m −1 . Previous studies reported that the lapse rate on snowfields and glaciers in the mid-latitude and polar region generally falls within the range during summer (Gardner et al., 2009;Wang et al., 2016). In Lr-exp, atmospheric conditions derived from WFDEI at each site were corrected using the method shown in Section 2.4 by applying different temperature lapse rates ( Figure S4 in Supporting Information S1). The initial cell concentration, growth rate, and carrying capacity were taken as 6.3 cells m −2 , 0.018 hr −1 , and the observed maximum cell concentration (cells m −2 ) for Rd-exp and Lr-exp, respectively. The experimental architecture for the sensitivity tests used in this study has been summarized in Table 3.
To evaluate snow algae model sensitivity to biological variables and atmospheric conditions, we used the two variables visible cell concentration of algal bloom (VCAB), and the reaching time to the algal bloom (RTAB). In this study, VCAB was defined as 5.0 × 10 5 cells m −2 because this was the minimum algal cell concentration for a red snow phenomenon observed on a Greenland glacier by Onuma et al. (2018) and were comparable to algal bloom concentrations observed in other areas (e.g., 2.2 × 10 6 cells m −2 at AT-A, 1.3 × 10 6 cells m −2 at TN-U, and 6.0 × 10 5 cells m −2 at PT-G). The time taken to reach the VCAB was defined as the RTAB. In the sensitivity tests shown in Table 3, if a calculated algal cell concentration did not reach the VCAB in a simulation, the date of the calculated maximum abundance was still defined as the RTAB.

Global Snow Algae Model Bio-MATSIRO and Experimental Design of Global Simulation
To evaluate snow algal seasonal growth changes and global distribution qualitatively, we incorporated the updated snow algae model (X MDR model) described in Section 2.1 into a scheme of snow physical processes in MAT-SIRO6. We named this version for snow algal simulation Bio-MATSIRO, using the same naming convention applied to the water isotope simulation, Iso-MATSIRO, by Yoshimura et al. (2006). Bio-MATSIRO can be used to calculate temporal changes in algal cell concentration (cells m −2 ) at regional and global scales, using atmospheric conditions near the land surface as input data. The model input data for simulations using Bio-MATSIRO in this work have been summarized in Table 2. In the global simulations shown in Section 3.8, we conducted a two-dimensional offline land simulation, using Bio-MATSIRO. For this, we conducted three global simulations, using Bio-MATSIRO with the WFDEI, GSWP3-FD, and CRUJRA reanalysis datasets. The horizontal resolution and calculation period for these simulations were 0.5° and from 1 January 1980 to 31 December 2014, respectively, to ensure that we had common horizontal resolutions and data periods for the three datasets. Land physical properties, such as snow water equivalent, were derived in advance from spin-up simulations (35 years total), using the same reanalysis dataset, and were used as initial conditions for the three simulations. The initial cell concentration and algal growth rate were the same as Ag-exp in this study (6.33 cells m −2 and 0.018 hr −1 , respectively) because these preliminary biometrics seemed to facilitate better model simulation performance. As there is little information available on carrying capacity, which is likely to vary between sites (snowfields or glaciers) and years, it was assumed to be 3.5 × 10 7 cells m −2 , which was the carrying capacity suggested by Onuma et al. (2018), for Greenlandic Glacier sites.

Simulation and Evaluation of Algal Cell Concentrations Achieved Using the Original Snow Algae Model
The algal growth simulations, Ag-exp, produced RTABs ranging from days 155-220; the date depended on their Northern Hemisphere locations, with mid-latitude site RTABs being generally earlier than those for polar sites. Specifically, the RTAB ranged from days 155-220 (early June to early August), at the Northern Hemisphere mid-latitude sites (N 30-60°; AT-A, EU-T, TN-U, NA-S, and JP-T), and from days 175-220 for the polar sites (N 60-90°; SV-S, SV-W, GL-Q, GL-A, AK-E, AK-K, and SB-S). The algal cell concentration, X M , at JP-T, which was calculated using snow melting duration only, showed no significant increase from days 1-70, and then increased, reaching VCAB on d 155 (RTAB = d 155; Figure 3). X M values for mid-latitude sites were found to reach the VCAB earlier than the polar sites. Similarly, X M at Southern Hemisphere mid-latitude sites (S 30-60°; NZ-A and PG-T) reached the VCAB earlier than it did at the single Southern Hemisphere polar site (S 60-90°; AP-L). The RTAB estimated from X MD , which was calculated from the duration of snow melting under daylight conditions, did not significantly differ from that estimated at the study sites using X M (Figure 3).
The simulation achieved using the original snow algae model showed algal cell concentration changes over time which agreed with in situ values for snowfields at GL-Q, AK-E, TN-U, and JP-T. For example, Ag-exp showed that the RTAB estimated for GL-Q from X M was days 215, which was the same date that the maximum cell concentration was observed at the site (Figure 3). The RTABs for sites AK-E, EU-T, TN-U, JP-T, and NZ-A also agreed well with the observations, whereas the changes in X MD over time did not significantly differ from those for X M , at any of the sites. The determination coefficients for temporal change in X MD (and X M ), compared with those for the observed algal cell concentrations, were 0.97 (P < 0.05), 1.00 (P < 0.05), 0.38 (P > 0.05), and 0.71 (P < 0.05), at GL-Q, AK-E, TN-U, and JP-T, respectively.

Effect of Daylight Length on Algal Bloom Simulation Results
The simulations showed that RTABs did not differ significantly in cases where daylight length was or was not taken into account, indicating that daylight length did not affect snow algal bloom timing to a significant extent. Ag-exp showed that GP M , which represents cumulative snow melting time, was longer than GP MD , excluding night time ( Figure S5 in Supporting Information S1). In the simulations, surface snow could melt at night when the daily surface air temperature was >10°C at mid-latitude sites ( Figure S1 in Supporting Information S1). However, even if the surface snow had melted, snow algae would be unlikely to grow during the night because they cannot photosynthesize without solar radiation. In experimental studies, the snow algae Cr. tughillensis and Cr. chenangoensis showed the greatest increased sexual reproduction under blue light with longest daylight scenario (light: dark 24: 0 hr; Hoham et al., 2000Hoham et al., , 2009, suggesting that snow algal growth depended on daylight length. Hence, GP MD , excluding nighttime, was expected to simulate algal growth better than GP M -however, there were no significant differences in the RTABs calculated using X M and X MD , at any site ( Figure 3). This indicated that snow did not melt during the night before algal blooming commenced. The reanalysis dataset showed that surface air temperatures often fell below freezing during the night in the early melting season, indicating that snow melting during the night would have been rare.

Simulation and Evaluation of Algal Cell Concentrations Achieved Using the Improved Snow Algae Model
RTABs simulated using the updated snow algae model ranged from days 180-240 in the Northern Hemisphere, showing that the snow algal bloom timing estimates were significantly later than those simulated using the original snow algae model, at all sites. For example, the algal cell concentration, X MDF , at JP-T, which was simulated using snow melting, daylight length, and interruption by new snow cover effects, kept X 0 from day 1-100, and then reached the VCAB on days 180 (RTAB = d 180; Figure 3). The RTAB ranged from day 180-255 (late June to early September) at the Northern Hemisphere mid-latitude sites, and from days 185-220 in the polar sites. In Southern Hemisphere, RTABs were 365, 355, and 40, at sites NZ-A, PG-T, and AP-L, respectively. These results showed that RTABs for X MDF was later than had been the case for X M and X MD . We also saw that there were no significant differences in the RTABs estimated using X MDF among the sites at different latitudes, unlike the case for those estimated using either X M or X MD .
Snow algal abundances simulated using the updated model produced seasonal change results which agreed better with observational data than those simulated using the original snow algae model. Ag-exp showed that X MDF at GL-Q started to increase on day 175 and then reached the VCAB on day 220. The temporal changes estimated using X MDF agreed with the site observation data better than those achieved using either X M or X MD (Figure 3). The determination coefficients for the temporal change in X MDF against the observation data were 0.97 (P < 0.05), 1.00 (P < 0.05), 0.52 (P > 0.05), and 0.79 (P < 0.05), for sites GL-A, AK-E, TN-U, and JP-T, respectively. These coefficients were slightly higher than those for the X M and X MD , whereas the RTABs calculated using X MDF agreed with the timing of the red snow phenomenon observed at the other sites, including SV-S, SV-W, EU-T, NA-S, NZ-A, PG-T, and AP-L.

Effect of Snowfalls on Algal Bloom Simulations
Updated model simulations which included the effects of snowfall and daylight length were more accurate than those produced using the original model, at most sites, suggesting that snowfall significantly affected algal growth, including algal bloom timing. The RTAB in the case of X MDF was later than it was the case of X MD , at all sites (Figure 3), whereas comparing temporal changes in the X MDF with those in the X MD showed that the rate of increase in X MDF declined due to frequent snowfall during the early melting season at the study sites. For example, daily surface air temperature at site JP-T first exceeded 0°C on day 90, but snowfalls continued to occur quite frequently, upto day 110. The algae gradually increased from day 90 to 110 at this site, in the case of X MD , whereas X MDF showed no increase during the same period. These results indicated that frequent snowfall events delayed the RTAB. Model simulations showed that increasing trends in X MDF agreed well with observed algal cell concentrations at sites GL-Q, AK-E, TN-U, and JP-T (Figure 3). Figure 6 indicates that the updated snow algae model was generally more accurate than the former model at many study sites, suggesting that incorporating snow cover on the algal growth surface had improved model accuracy. Although the effect of snowfall on an algal surface had been previously reported (Onuma et al., 2016;Tanaka, 2016), it had not been quantified in field observations. Further periodical observations during the early snow melting season are now needed to further improve spring to summer snow algal bloom prediction accuracy.

Sensitivity of Algal Growth Simulation to Initial Cell Concentration and Growth Rate
Testing X 0 -exp sensitivity to the initial cell concentration showed that the RTAB was approximate 15-30 days earlier where the initial cell concentration was 100-fold greater than the original concentration (1.0 cells m −2 , Figure 4a). For example, the RTABs in the simulation of the minimum initial cell concentration (1.0 cells m −2 ) were on days 215 and 195, at AK-E and JP-T, respectively. Simulation results showed that the AK-E RTABs were on days 200, 180, and 175, for the cases in which the minimum initial cell concentrations increased by a factor of 10, 100, and 1,000, respectively. Under similar circumstances, the RTABs for JP-T were on days 180, 160, and 145. X MDF did not reach the VCAB before the disappearance of snow in either GL-A or AP-L, in any of the simulations.
The test of μ-exp sensitivity to algal growth rate showed that, in the case of a 10% greater growth rate, the RTAB was approximately 10 days earlier than the original rate (Figure 4b). RTABs simulated using a growth rate of 0.018 hr −1 (original case) were 210 and 180 days at AK-E and JP-T, respectively. The RTABs at AK-E were simulated to be days 200, 195, and 190, for growth rates 10%, 20%, and 30% greater than that of the original rate, respectively. Under similar circumstances, the JP-T RTABs were simulated to be days 170, 160, and 155. X MDF did not reach the level of algal blooming at GL-A, SB-S, EU-T, NZ-A, or AP-L, in any of the simulations.
Sensitivity testing suggested that the RTAB range was likely to be related to site latitude. For example, the RTAB at JP-T (36°N) ranged from days 145-195, whereas, at SV-W (77°N), it ranged from days 180-200 in X 0 -exp ( Figure 4a). Similarly, the RTAB at JP-T ranged from days 160-250, whereas it ranged from days 180-200 in μ-exp ( Figure 4b). These differences in the RTAB among the sites could be explained by their different GP MDF results. GP MDF accumulated during the daytime only at JP-T, whereas at site SV-W, it also accumulated during the night ( Figure S5 in Supporting Information S1). The sensitivity tests showed similar trends at the other sites. At polar sites, snow algae are likely to grow for a shorter period because of the longer daylight length during summer, leading to a shorter range for the RTAB in each sensitivity test. For example, the RTAB simulated with the initial cell concentration of 694 cells m −2 was 15 days earlier than that simulated using 6.3 cells m −2 at SV-W and 30 days earlier in the case of JP-T. In μ-exp, the RTAB with a growth rate of 0.018 hr −1 was 10 days earlier than that simulated with an initial cell concentration of 0.016 hr −1 , at SV-W, whereas it was 25 days earlier in the case of JP-T. Although differences in biological parameters could lead to uncertainty in algal bloom timing, the uncertainty may be smaller in polar regions than in mid-latitude regions. The results here suggested that the timing estimates for the red snow phenomenon achieved using the updated snow algae model were more reliable at polar sites.
Testing biological parameter sensitivity suggested that the RTAB estimates in each test significantly change within the range of initial cell concentrations and growth rates reported by the previous study. In this study, we conducted simulations using the updated snow algae model with an initial cell concentration of 6.33 cells m −2 , and a growth rate of 0.018 hr −1 . However, these parameters may have introduced some uncertainty, as they were derived using in situ Greenland glacier data (Onuma et al., 2018). S. nivaloides blooming was previously reported as appearing on days 180 and 215 at two different sites, with the initial cell concentration and the growth rate at the former site (694 cells m −2 and 0.018 hr −1 ) higher than those at the latter site (6.33 cells m −2 and 0.016 hr −1 ). Differences in biological parameters would cause uncertainty in the timing of snow algal blooming, and thus, we quantified the range of RTAB uncertainty which could be attributed to the biological parameters observed in the previous study. The range was obtained from simulation results (X 0 -exp and μ-exp) using the lowest (6.33 cells m −2 and 0.016 hr −1 ) and highest (694 cells m −2 and 0.018 hr −1 ) reported biological parameters ( Figure 5).
RTABs simulated using the lowest and highest parameters are shown in Figures 5a and 5b, respectively, indicating that the RTABs in Figure 5b were 10-30 days earlier than those in Figure 5a. Snow algal blooming was generally observed on the day between the RTABs in Figures 5a and 5b, suggesting that the RTAB simulated with the updated snow algae model changed by approximately 10 days, compared with the observed blooming date of the snow algae. Although snow algae initial cell concentrations and growth rates may depend on the amount of mineral dust supplied from the atmosphere and on nutrient conditions (nitrogen and phosphorus), as suggested by Onuma et al. (2016Onuma et al. ( , 2018, the major factors affecting the biological parameters remain uncertain. Field data on snow algae biological parameters from various snowfields and glaciers worldwide are needed. However, our  simulations suggested that the updated snow algae model could forecast snow algal bloom timing with an accuracy of approximately 10 days.

Algal Bloom Simulation Uncertainties Caused by Different Atmospheric Conditions
Testing Rd-exp sensitivity to atmospheric conditions showed that the RTAB estimate results significantly varied among the reanalysis datasets of the global atmosphere, even at the same site. The minimum, median, and maximum RTABs, derived from simulations using the WFDEI, GSWP3-FD, and CRUJRA datasets, can be seen for each site in Figure 6a-which also shows the reanalysis dataset uncertainties for the RTABs. The RTABs varied significantly among the datasets at sites SV-S, GL-A, SB-S, AT-A, EU-T, TN-U, and NA-S. The differences in the RTABs estimated between the datasets were equal to 100 days at EU-T and just 5 days at site JP-T.
Rd-exp showed that sensitivity to meteorological conditions was higher in the Asian high mountain areas than in polar snowfields, suggesting that RTAB estimation uncertainty was larger in areas where precipitation mainly occurred in summer. RTABs simulated with the updated snow algae model using the WFDEI, GSWP3-FD, and CRUJRA reanalysis datasets varied greatly at some sites ( Figure 6a). Notably, the difference in the RTAB between the 25th and 75th percentiles was 50 days at sites EU-T and TN-U. The differences were smaller, approximately 20 days, at some alpine sites (SB-S, AT-A, and NA-S) and polar maritime snowfield sites (SV-S, SV-W, and GL-A). This RTAB difference was probably due to the frequency of snowfall during summer in each dataset. As reported previously, most Asian high mountain glaciers are characterized by summer accumulation due to the influence of the Asian monsoon (Fujita & Ageta, 2000;Fujita, 2008;Sakai & Fujita, 2017). According to meteorological conditions derived from the reanalysis datasets, study sites SB-S, AT-A, EU-T, and TN-U could probably be classified as summer accumulation-type glaciers ( Figures S1 and S3 in Supporting Information S1). The higher sensitivity of the RTAB to atmospheric conditions was probably caused by frequent summer snowfalls at these sites, which would greatly affect the algal growth in the X MDF . Furthermore, a previous study showed that precipitation amounts derived from the reanalysis dataset of the global atmosphere still have levels of uncertainty in high elevation and polar maritime areas (Weedon et al., 2014). Large RTAB variations may also be influenced by dataset accuracy at specific sites-although there were no significant differences in the RTABs estimated using the different datasets at sites GL-Q, AK-E, and JP-T. The determinant coefficient in the algal cell concentrations achieved by simulation and observation at these sites was >0.79 (Figure 3), suggesting that the updated snow algae model was capable of reproducing red snow bloom timing in these regions, even if the dataset contained uncertainties regarding atmospheric conditions. Notably, the updated snow algae model performed very reasonably in reproducing snow algal bloom timing at sites where the reanalysis data of the global atmosphere were highly accurate.
Testing Lr-exp sensitivity to temperature lapse rate showed that there was no significant difference in the RTAB estimate results among the different lapse rates, even at the same site. The minimum, median, and maximum RTABs, derived from simulations using the WFDEI with the ranging from −8.0 × 10 −3 to −5.0 × 10 −3 K m −1 , can be seen for each site in Figure 6b-which also shows the temperature lapse rate uncertainties for the RTABs. The RTABs varied significantly among the lapse rates at site GL-Q only. The differences in the RTABs estimated between the lapse rates were within 10 days at the other sites.
Lr-exp showed that sensitivity to temperature lapse rates was significantly higher in the polar high elevation site, suggesting that RTAB estimation uncertainty was larger in areas where the surface air temperature was around the melting point throughout summer. RTABs simulated with the updated snow algae model using the WFDEI significantly varied among the different lapse rates at GL-Q only (Figure 6b). The difference in the RTAB between the 25th and 75th percentiles was 20 days at GL-Q, whereas it was below 10 days at the other sites. There was no significant difference in the surface air temperature after elevation correction among the different lapse rates at most sites ( Figure S4 in Supporting Information S1). This is probably due to that the elevation information of the WFDEI grids near the validation sites was similar to that of the validation sites. Although the surface air temperature in GL-Q did not significantly vary among the lapse rate, the air temperature was around 0°C throughout summer. In such atmospheric conditions, the uncertainty of the algal growth simulation would be higher.
The model needs to be improved if it is to reproduce red snow blooming accurately on a global basis-especially in summer accumulation-type glaciers. At such sites, model validation using in situ observational meteorological conditions would be necessary. It has also become apparent that the level of snow algal bloom sensitivity to meteorological conditions in summer may vary depending on seasonal precipitation patterns. Further, the snow algae model could be a useful tool for revealing snow algal growth sensitivity to meteorological conditions.

Other Possible Factors Affecting Red Snow Algal Blooming
The current status of snow algal growth numerical modeling in snowfields has been summarized in this section, with aspects of the updated snow algae model which still require improvement being identified. The updated model can estimate algal bloom timing better than what had been achievable in simulations using the previous model (Onuma et al., 2018) at comparable sites ( Figure 5). As noted in Section 3.4, this was probably due to incorporating the effect of snowfall on the algal abundance in surface snow into the model. The result also showed, however, that algal cell concentrations simulated using the updated model underestimated observed concentrations at the study sites, resulting in algal blooming not appearing during the thaw season at the polar sites GL-A, SB-S, and AP-L. This may be attributed to uncertainties in both the biological parameters and the reanalysis data (especially those for snowfall amount and frequency) at the sites. Other biological process issues may also be contributing to these underestimations. Previous studies have suggested that snow algae motile cells could swim up to the snow surface from the soil, or from the ice below the surface snow (Müller et al., 2001;Remias, 2012).
The snow depth at site AP-L ranged from 5-20 cm during the observation period, and thus, as Müller et al. (2001) suggested that such motile cell vertical movement occurred at snow depths <40 cm, snow algal cells originating from the ground may have contributed to underestimating the simulated algal cell concentration at AP-L. In addition to the supply of the algal cells from the soil, previous studies reported that algal cells are supplied on the snow surface from the atmosphere (Onuma et al., 2016(Onuma et al., , 2018. In this study, the snow algae model neglects the effect of these supplied algal cells on the simulated algal abundance. Although the snow algae model could simulate temporal changes in cyst cell concentrations, some of these cells might transform into vegetative cells during the thaw (Remias et al., 2010). Moreover, the algal growth rate and the carrying capacity may vary among the snow algal species. Further field observation and cultivation are needed to quantify such algal stage changes so that they can be incorporated into the snow algae model. Because the model equations and biological parameters are simple in this study, the latest information regarding snow algal growth should be incorporated into the snow algae model in future.

Global Simulation of Red Snow Algae Blooms
The global simulation established using Bio-MATSIRO showed that snow algae grew from spring to summer in both hemispheres and that their blooming sites were generally consistent with red snow sites reported previously, suggesting that Bio-MATSIRO had the potential to reconstruct snow algal blooms at the global scale. The global distribution of snow algal blooming was derived from the monthly means of the X MDF (1980-2014 climatological mean), which were the algal cell concentrations simulated with Bio-MATSIRO using each reanalysis dataset. Red snow blooming (X MDF , reanalysis dataset of the global atmosphere: WFDEI) distribution for each month's snow cover can be seen in Figure 7. In the Northern Hemisphere, snow algae gradually increased from March to August, and their blooming area gradually extended from mid-to high-latitude areas. The area covered by the red blooms then extended southward because snowfall events frequently occurred in high latitudes. Red snow algal blooming simulated with Bio-MATSIRO appeared between June and August, and the blooming area agreed well with the red circles (red snow phenomenon reported by previous studies) shown in Figure 1c. Red snow in European and North American snowfields had been reported previously over periods extending from summer to early autumn (Remias et al., 2016;Thomas & Duval, 1995), and the reported areas were consistent with simulated red snow algal bloom areas. In the Southern Hemisphere, the model simulation presented red snow algae gradually increasing from September to March, especially in Patagonia and on the Antarctic Peninsula. The red snow phenomenon has been reported as occurring between November and March in these regions (Gray et al., 2020;Takeuchi & Kohshima, 2004). Interestingly, the red snow algal bloom simulation hardly showed any occurrences in either N Russia or NW Canada-and there have actually been few reports of this phenomenon there. Frequent snowfall events during the snow melting season would interrupt snow algal bloom appearance in these regions because X MDF drastically decreased to X 0 in August. Simulations using GSWP3-FD and CRUJRA showed spatial and seasonal changes in X MDF similar to those simulated using WFDEI ( Figures S6 and S7 in Supporting Information S1), with these results suggesting that Bio-MATSIRO had the potential to reproduce seasonal and geographical changes in red snow algal abundance globally.
Because there is still uncertainty of the biological parameters (initial cell concentration X 0 , growth rate μ, and carrying capacity K) for each region used in this study, we additionally conducted sensitivity tests to biological parameters at the global scale. Figure 8 shows RTAB simulated using Bio-MATSIRO with various biological parameters. The sensitivity ranges of the X 0 and μ were based on the field observation by Onuma et al. (2018). The range of K was assumed from the maximum algal cell concentration shown in Table 1. Although the result indicates that grids, where red snow appeared, were expanded in the case of using greatest X 0 , μ, and K, there was no significant difference in the expanded area among the parameters. However, the RTAB was generally 15-30 days earlier where the X 0 was 100-fold greater than the least concentration (6.33 cells m −2 , Figure 8c). On the other hand, the RTAB was only 5-10 days earlier in the cases of greatest μ and K. These results suggest that the sensitivity to the initial cell concentration is higher rather than the other parameters. Onuma et al. (2018) suggested that the initial cell concentration depends on the mineral dust abundance supplied from the atmosphere. Quantification of the supplied mineral dust abundance in the surface snow is likely to be a key point to quantify the initial cell concentration using a numerical model. Another study suggests that red snow blooming comprise mainly cosmopolitan phylotypes but also includes endemic organisms (Segawa et al., 2018). Such a process may be also important to reveal their growth and initial cell concentration worldwide. In addition to the initial cell concentration, further observation is necessary to reveal a major factor affecting the algal growth rate and the carrying capacity. As Figure 8 showed the uncertainty of the biological parameters, the algal cell concentrations derived from our simulations contained some uncertainty levels in this study. The uncertainty may lead to that the model simulation estimates the timing of the red snow appearance earlier or later than the real world. Further observational and modeling studies are necessary to improve the snow algae model. Periodical observation is important to quantify supplied algal cells from the atmosphere or soil, nutrient conditions in the snowpack, transformation rate into vegetation cells. Satellite observations of red and green snow algae blooms have been conducted recently (Ganey et al., 2017;Gray et al., 2020;Hisakawa et al., 2015;Huovinen et al., 2018;Khan et al., 2021), and validating model simulations with satellite observations would also be useful, at the glacial or regional scales. In addition, global simulation of red snow algal blooming, using a land surface model, could provide an important contribution to understanding climate change effects on snow and ice distribution in time and space.

Conclusions
We updated the existing snow algae model based on observational data from 15 snowfields and incorporated it into a land surface model to quantify time and space changes in snow algal abundance worldwide. The existing snow algae model (Onuma et al., 2018) could simulate temporal changes in the abundance of red snow algae in surface snow, using snow temperature, but up until now had only been applied to simulate the abundance of Greenlandic glacier snowpack. In this study, the effects of daylight length and snowfall rate on algal cell abundance were incorporated into the model, and the revised model simulations achieved good agreement with observations at snowfields worldwide, from polar to mid-latitude areas-particularly in regions with fewer summer snowfalls. Based on these encouraging results, we incorporated the updated snow algae model into a land surface model and conducted a global snow algal simulation, using Bio-MATSIRO. This simulation produced results showing prominent algal blooms taking place in areas generally consistent with regions where the red snow phenomenon had been reported in either in situ or satellite observations. Our simulations suggested that Bio-MATSIRO has the potential to simulate temporal and spatial changes in red snow algal abundance and to predict the timing and coverage of the red snow phenomenon from the past to the future. The simulations also suggested that the timing of red snow appearance is likely to be the balance between the duration of the snow melting and the timing of the snowfall event, meaning that the habitat of snow microbes is closely linked to the climatic environment.
Snow algal distribution may be a key to revealing the geographic specifications of microbes worldwide (Lutz et al., 2016;Procházková et al., 2019;Segawa et al., 2018;Zawierucha & Shain, 2019). The snow algae model may be useful not only for providing such biological information but also for quantifying snow algal contributions to thaw events worldwide (by reducing snow surface albedo), revealing trends over time. This contribution has been reported in many studies as a bioalbedo effect and has been quantified using in situ observations, satellite observations, and numerical simulations (Aoki et al., 2013;Cook et al., 2017;Ganey et al., 2017;Gray et al., 2020;Lutz et al., 2016;Mauro et al., 2017;Onuma et al., 2020;Painter et al., 2001;Takeuchi, Dial, et al., 2006;Thomas & Duval, 1995). Because these studies focused on specific mountains or glaciers, comprehensive bioalbedo effects should be investigated at the global scale in the future using numerical simulation. To simulate bioalbedo effects at such a scale, a numerical bioalbedo model that has the capacity to calculate snow albedo (with the effect of snow algae included) needs to be incorporated into land surface models and climate models. The snow algae model coupled with the bioalbedo model might enable to project of the effect of snow algal blooming on climate through variations in snow albedo. Although further observations and simulations would further improve snow algae models, our study has provided an important first step toward revealing the global geographic characteristics of snow algae and their contribution to snow melting.