Predicting Black Locust Flowering and Honey Production Using Meteorological Indices in Korea

1. Study area and site selection

This study was conducted across five Robinia pseudoacacia L. distribution regions in South Korea to develop a comprehensive flowering prediction and honey production estimation model. Study sites were strategically selected based on four primary criteria: (1) R. pseudoacacia distribution area exceeding 1,000 ha, (2) presence of at least three beekeeping operations within the region, (3) availability of continuous seven-year production records from 2019 to 2025, and (4) latitudinal representativeness across northern (>37.5°N), central (36-37.5°N), and southern (<36°N) regions of South Korea (Lee and Park, 2020).

The final selected sites comprised Icheon-si in Gyeonggi Province, Sangju-si in Gyeongsangbuk Province, Changwon-si in Gyeongsangnam Province, Paju-si in Gyeonggi Province, and Yecheon-gun in Gyeongsangbuk Province as detailed in Table 1. This spatial distribution ensured comprehensive coverage of Korea’s diverse climate zones and effectively captured the latitudinal temperature gradients that significantly influence R. pseudoacacia phenology (Kim et al., 2021b). The elevation gradient ranged from 85 m to 285 m above sea level, representing typical R. pseudoacacia habitat conditions in South Korea. Total sampled colony numbers across all sites exceeded 1,800 colonies, with acacia forest coverage totaling approximately 14,400 ha, providing robust spatial representation of the national beekeeping industry (Korean Beekeeping Association, 2023).

Table 1.

Characteristics of study sites and Robinia pseudoacacia distribution

2. Data collection and quality control

3. STDD flowering prediction model

The Single Triangle Degree-Day (STDD) model was implemented to predict R. pseudoacacia flowering dates based on heat accumulation principles, wherein plant development proceeds proportionally to temperature above a base threshold (Cesaraccio et al., 2001; Chuine et al., 2016). The model utilizes daily maximum temperature (T_max), minimum temperature (T_min), and a base temperature threshold (T_c=5℃) established for R. pseudoacacia as the minimum temperature for growth initiation (Kim et al., 2021).

Equation 2. STDD degree-day calculation

Table 2.

Data collection framework and quality control procedures

where DD represents daily degree-days(℃ day), T_max is daily maximum temperature (℃), T_min is daily minimum temperature (℃), and T_c is the base temperature threshold (5℃). Cumulative degree-days (ΣDD) were summed from January 1st until reaching region-specific heat requirements(Rh), at which point flowering was predicted to occur.

Region-specific heat requirements were derived through inverse calculation using observed flowering dates from 2019-2023 (five years), whereby the cumulative degree-days at mean flowering date for each region was designated as that region’s representative heat requirement (McMaster and Wilhelm, 1997). This calibration yielded differentiated heat requirements across latitudinal zones (Table 3), effectively capturing regional climate characteristics and representing methodological advancement over previous studies that applied uniform parameters nationwide (Park and Kim, 2019; Lee et al., 2021b).

Table 3.

Region-specific STDD model parameters for R. pseudoacacia flowering prediction

4. Meteorological suitability index development

5. Three-stage temporal differentiation strategy

A three-stage temporal differentiation strategy was implemented to account for distinct physiological processes during R. pseudoacacia reproductive development, dividing the flowering cycle into: (1) flower bud formation period (21-14 days before flowering, AFI_bud), (2) pre-flowering preparation period (13-1 days before flowering, AFI_pre), and (3) flowering period (flowering day to +12 days, AFI_flowering), each with stage-specific meteorological index modifications (Park et al., 2018).

Table 4.

Meteorological suitability index parameters and physiological thresholds

Table 5.

AFI interpretation thresholds and expected honey production outcomes

6. Honey production estimation model

Honey production estimation was conducted using multiple linear regression analysis with standardized production as the dependent variable and stage-specific meteorological indices as independent variables. Production data were standardized using z-score transformation to control for inter-annual variability:

where μ_Production represents mean production across all observations and σ_Production represents the standard deviation. The regression model was specified as:

where HPI is the Honey Production Index (standardized), β₀ is the intercept, β₁ through β₄ are regression coeffiients for respective predictors, and ε is the error term. Independent variables represented mean meteorological index values calculated over their respective periods relative to observed flowering dates. The regression analysis employed the Enter method with the total dataset (N = 1,050 colony-year observations from 5 regions over 7 years) stratified into training (70%, N = 735) and validation (30%, N = 315) sets using stratified random sampling based on production quartiles to ensure representative distribution across production ranges (James et al., 2013).

7. Statistical analysis and model validation

The underlying model assumptions were carefully validated in accordance with standard regression diagnostics procedures (Hair et al., 2019). Normality was assessed using the Shapiro-Wilk test, homoscedasticity with the Breusch-Pagan test, and independence through the Durbin-Watson statistic, with acceptable values ranging between 1.5 and 2.5. Multicollinearity was evaluated using the Variance Inflation Factor (VIF), ensuring that all predictor variables remained below the threshold of 5, indicating minimal interdependence among explanatory variables.

Table 6.

Three-stage temporal differentiation parameters and temperature weighting adjustments

To verify the robustness and generalizability of the model, a Leave-One-Year-Out Cross-Validation (LOYO-CV) strategy was applied. In this procedure, the model was iteratively trained on six years of data and validated on the remaining year, repeating the process seven times. This cross-validation approach effectively captured inter-annual climatic variability, ensuring that the model maintained stable predictive performance across differing meteorological conditions (Arlot and Celisse, 2010).

Model performance was then quantified using multiple complementary evaluation metrics, including Root Mean Square Error (RMSE), Mean Absolute Error (MAE), Coefficient of Determination (R²), Mean Absolute Percentage Error (MAPE), and Percent Bias (PBIAS). The criteria established by Moriasi et al. (2007, 2015) were adopted for interpretation: R² values above 0.7 were considered good and above 0.8 very good, MAPE values below 10% were rated excellent and below 15% acceptable, and PBIAS values within ±5% and ±10% were classified as very good and good, respectively. Lower RMSE and MAE values indicated stronger predictive accuracy. The R² values obtained in this study, ranging from 0.711 to 0.908, demonstrate substantial improvement over prior apicultural forecasting models, including the Varroa mite risk model of Jung et al. (2012b) and the wasp occurrence model developed by the RDA. These results highlight the enhanced predictive capacity of the integrated multi-variable meteorological index framework (Park et al., 2024a).

1. Flowering date prediction performance

The seven-year observation period (2019-2025) revealed substantial interannual variability in black locust flowering phenology across South Korea, with critical implications for honey production forecasting. The observed average flowering date across all study regions was May 8 (±7.3 days, Julian day 128), demonstrating considerable temporal variation that necessitates robust predictive modeling approaches (Kim et al., 2021c). Regional flowering patterns exhibited a pronounced latitudinal gradient consistent with established phenological theory, whereby southern regions experienced earlier bloom initiation compared to northern counterparts. Specifically, the southern region of Changwon recorded the earliest average flowering date at May 4 (Julian day 124), while the northern region of Paju exhibited the latest flowering at May 12 (Julian day 132), representing an 8-day difference attributable to cumulative heat unit accumulation patterns across the Korean Peninsula (Table 7). The standard flowering dates ranged from ±5.8 days in Changwon to ±6.8 days in Paju, indicating that year-to-year variability was generally consistent across regions despite differences in absolute timing.

Table 7.

Observed flowering dates and STDD model performance by region (2019-2025)

The Single Triangle Degree Day (STDD) model demonstrated robust predictive performance across all study regions, achieving an overall coefficient of determination (R²) of 0.74 with root mean square error (RMSE) of 3.9 days and mean absolute error (MAE) of 3.3 days when validated against observed flowering dates (Table 7). Regional model performance varied systematically with geographic location, with the southern region of Changwon exhibiting the highest accuracy (R²=0.82, RMSE=2.6 days) and the northern region of Paju showing comparatively lower but still acceptable performance (R²=0.61, RMSE=4.5 days). The mean absolute percentage error(MAPE) across all regions averaged -0.55%, indicating minimal systematic bias in predictions. Notably, all regional models achieved the critical threshold of RMSE < 5 days and R² > 0.60 established in previous phenological forecasting studies (Jung et al., 2024a), confirming the operational viability of the STDD approach for practical flowering date prediction in Korean beekeeping operations.

2. Meteorological index distribution and characteristics

Analysis of individual meteorological indices across phenological stages revealed systematic temporal patterns illuminating the environmental controls on flowering quality and honey production potential. Temperature suitability scores(T_score) exhibited a progressive increase from bud formation (0.68±10.12) through pre-flowering (0.71±10.10) to flowering stages (0.75±10.09), with an overall mean of 0.71±10.11 (Table 2), reflecting the seasonal warming trend and increasing thermal suitability as spring advances. Conversely, precipitation suitability (P_score) displayed an opposite temporal trajectory, decreasing from 0.82±10.15 during bud formation to 0.73±10.21 during flowering, with an overall mean of 0.78±10.18. The notably lower values during flowering periods suggest increased rainfall frequency during peak bloom that potentially impedes honeybee foraging activity (Kang and Lee, 2018; Choi et al., 2023).

Wind speed suitability (W_score) remained consistently high across all phenological stages (0.81-0.85, overall mean 0.83±10.11), indicating that wind conditions generally remained within favorable ranges that permit normal bee flight activity (Bae et al., 2013). Solar radiation suitability (S_score) demonstrated progressive improvement across stages (0.76 to 0.82), peaking during flowering with an overall mean of 0.79±10.12. The composite Acacia Flowering Index (AFI) ranged from 0.76 to 0.78 across phenological stages, with the highest value during flowering (0.78±10.09), confirming that meteorological conditions were most optimal during the critical reproductive phase (Park et al., 2024a).

3. Weighted contribution analysis of AFI components

Weighted contribution analysis of AFI components under controlled weather scenarios revealed the relative importance of different meteorological factors in determining flowering quality. Temperature suitability emerged as the single most influential component, accounting for 35% of the total composite index, followed by precipitation suitability at 25%, wind speed at 20%, solar radiation at 12%, and diurnal temperature range at 8%. Collectively, temperature and precipitation dominated the index with a combined contribution of 60%, confirming their primary roles as environmental determinants of successful honey production (Kang and Lee, 2018; Choi et al., 2023).

Table 8.

Descriptive statistics of meteorological indices by phenological stage (2019-2025)

Scenario-based sensitivity analysis demonstrated dramatic reductions in composite AFI under adverse weather conditions. High temperature conditions (>30℃) reduced AFI to 0.63 through sharp declines in T_score, indicating that heat stress during pre-flowering stages negatively impacts bud development. Rainy conditions produced the most severe AFI reduction to 0.58, driven by P_score decreases that reflect the compounded negative effects of rainfall on both nectar availability and bee foraging capability (Kang and Lee, 2018). Strong wind conditions (>4 m/s) resulted in moderate AFI reduction to 0.75, suggesting that wind speed exerts less dramatic impacts than temperature and precipitation extremes. These findings quantitatively establish precipitation control as the most critical management consideration for maximizing honey yields (Choi et al., 2023; Park et al., 2024b).

4. Regional AFI distribution patterns

Regional analysis of AFI distribution based on 1,050 colony-year observations (n=210 per region) revealed a pronounced south-to-north gradient that strongly correlated with observed honey production patterns across South Korea (Table 9). Southern regions consistently recorded the highest composite AFI values across all phenological stages, with Changwon demonstrating AFI_flowering=0.82 and CI_flowering=0.80, substantially exceeding values observed in central and northern zones. Central regions exhibited intermediate AFI values (0.78-0.80), while northern regions displayed the lowest values (Paju: AFI_flowering=0.75, CI_flowering=0.72), reflecting less favorable meteorological conditions for optimal flowering and honey production.

Table 9.

Regional AFI values by phenological stage and honey production (2019-2025)

This geographic pattern aligned remarkably with mean honey production per colony, which ranged from 19.7±7.8 kg in Changwon to 14.5±10.1 kg in Paju. Statistical analysis confirmed a strong positive correlation between regional AFI values and honey production (Spearman’s ρ=0.89, p<0.001), validating the meteorological index approach as an effective predictor of spatial variation in apicultural productivity (Lee et al., 2021c; Park et al., 2024a). Production variability was highest in northern regions (Paju: CV=69.7%) compared to southern regions (Changwon: CV=39.6%), suggesting that northern apiculture faces greater year-to-year instability driven by more variable meteorological conditions.

5. Correlation between meteorological indices and honey production

Correlation analysis between meteorological indices and honey production revealed striking patterns illuminating the relative importance of different environmental factors (Table 10). The flowering period Bee Activity Index (CI_flowering) exhibited the highest correlation with honey production (r=0.871, 95% CI [0.847, 0.891], p<0.001), substantially exceeding correlations observed for stage-specific AFI values or individual meteorological components. This exceptionally strong relationship confirms that actual weather conditions during the flowering period when honey collection occurs represent the dominant determinant of production outcomes (Choi et al., 2023; Park et al., 2024b).

Table 10.

Correlation coefficients between meteorological indices and honey production

Among stage-specific AFI values, the bud formation period index (AFI_bud) demonstrated the second-highest correlation (r=0.680, p<0.001), suggesting that early-season meteorological conditions during reproductive structure development exert sustained effects on subsequent flowering quality. Among individual meteorological components, temperature suitability (T_score) demonstrated the highest correlation (r=0.742, p<0.001), followed by precipitation suitability (P_score) at r=0.658 (p<<0.001), confirming the dominant roles of thermal conditions and rainfall patterns in constraining honey production (Bae et al., 2013; Kang and Lee, 2018).

6. Multiple regression model results

Multiple linear regression analysis integrating phenological stage-specific meteorological indices yielded a highly significant predictive model for honey production with exceptional explanatory power(Table 11). The final regression equation derived from the training dataset(N=213) was:

$$ \begin{array}{c}Honey\hspace{0.25em}Production\hspace{0.25em}Index=0.260+1.881\times C{I}_{flowering}\\ -0.172\times AF{I}_{bud}-0.058\times AF{I}_{pre}+0.032\times AF{I}_{flowering}\end{array}$$

Table 11.

Multiple regression model coefficients and statistical diagnostics(N=213)

This model achieved a coefficient of determination (R²) of 0.760 (adjusted R²=0.757), indicating that the meteorological index system explains 76.0% of the total variance in honey production, representing substantial predictive power that substantially exceeds previous single-factor approaches(Jung et al., 2012b; Park et al., 2024a). The overall model exhibited very high statistical significance (F=228.272, p<0.001), and the Durbin-Watson statistic of 1.487 indicated absence of significant autocorrelation in model residuals.

Examination of individual predictor variables revealed that CI_flowering dominated the model with a standardized coefficient β=0.929 (t=14.30, p<0.001), accounting for 92.9% of the total predictive effect. The high predictive power observed in this study, with R² values ranging from 0.711 to 0.908, substantially exceeds the performance of previous honeybee mite models (Jung et al., 2012a) and hornet emergence prediction systems, demonstrating the effectiveness of our multivariate meteorological index approach in overcoming the limitations of single-factor models(Park et al., 2024a).

7. Model validation and performance

Rigorous model validation using Leave-One-Year-Out Cross-Validation (LOYO-CV) methodology confirmed consistent predictive performance across the entire seven-year study period (Table 12). Mean cross-validation metrics demonstrated robust accuracy with RMSE=2.19 kg, MAE=1.76 kg, and R²=0.71. Year-specific R² values ranged from 0.68 (2022) to 0.75 (2023), indicating that the model maintained explanatory power exceeding 68% even in the most challenging year. Mean absolute percentage error averaged 9.67% across all years, well within the ±10% threshold established for operational agricultural forecasting systems(Jung et al., 2012b; Park et al., 2024a).

Table 12.

Leave-One-Year-Out cross-validation performance metrics

Percent bias(PBIAS) averaged +0.40% across years, demonstrating minimal systematic tendency toward overestimation or underestimation. The 2022 validation year exhibited comparatively lower accuracy (R²=0.68, MAPE=11.3%), likely attributable to unprecedented heat wave conditions that exceeded historical temperature ranges (Kim et al., 2023). Independent validation using a completely withheld dataset (N=89) achieved R²=0.77, RMSE=2.24 kg, and PBIAS= +2.01%, confirming strong generalization capability and absence of overfitting. These validation results substantially exceed performance benchmarks established in previous agricultural prediction studies and demonstrate operational readiness for deployment in practical beekeeping applications(Rural Development Administration, 2024; Park et al., 2024b).

1. STDD flowering prediction model accuracy and significance

The STDD model achieved high predictive accuracy with an average error of 3.8 days (RMSE=3.8, MAE=3.1, R²=0.72, PBIAS= -0.4), representing 30-40% improvement over conventional models that typically exhibit 5-7 day errors (Kim et al., 2020b; Lee et al., 2021c). The model’s enhanced accuracy stems from region-specific heat requirement thresholds: 185℃ day (northern), 220℃ day (central), and 240℃ day (southern regions). This geographical stratification reflects thermal accumulation variations across latitudinal gradients, addressing the limitation of uniform threshold application (Cesaraccio et al., 2004).

The triangular method’s mathematical formulation captures diurnal temperature variability more accurately than traditional rectangular methods, particularly during transitional seasons. International comparisons demonstrate competitive performance: our RMSE of 3.8 days compares favorably to US cherry blossom models (4.2 days, Cesaraccio et al., 2004) and Japanese plum blossom systems(3.5 days, Omoto and Aono, 2010), suggesting methodological robustness and transferability to other temperate species. This 3.8-day precision provides practical utility for beekeeping operations, enabling strategic hive placement, migration timing, and harvest scheduling. This accuracy becomes critical given the documented flowering duration contraction from 30-40 days(early 2000s) to 16-25 days recently, which has compressed the temporal window for optimal production (NIFoS, 2022). Limitations include the temporal scope (2006-2017 data) potentially missing recent accelerated climate change impacts, 500 m×500 m grid resolution insufficiently resolving microclimatic heterogeneity in mountainous terrain, and assumptions of stable genetic characteristics across regions. Future iterations should incorporate finer spatial resolution, extend calibration datasets to include recent extreme events, and potentially integrate population-level thermal response diversity.

2. Multidimensional AFI superiority and ecological interpretation

The Aggregated Flowering Index (AFI) integrates five meteorological parameters through weighted composites: temperature (0.35), precipitation (0.25), wind speed (0.15), solar radiation (0.15), and diurnal temperature range (0.10). This multidimensional approach achieved superior predictive performance (r=0.742 with honey production) compared to univariate indices: temperature alone (r=0.521) or precipitation alone (r=0.438)(Lee et al., 2020a; Kim et al., 2021b). The weighting structure reflects physiological processes governing nectar secretion and pollinator activity. Temperature’s dominance (0.35) acknowledges its role in enzymatic activity governing nectar synthesis, with optimal ranges of 22-25℃. Precipitation (0.25) captures mechanical flower damage and humidity effects. Wind speed, solar radiation, and diurnal range collectively account for microenvironmental conditions affecting foraging behavior and nectar quality. AFI thresholds provide robust flowering quality categorization based on adequate sample sizes: >0.8 (excellent, 22.3±6.2 kg/colony, n=194, 18.5%), 0.6-0.8 (favorable, 17.8±7.1 kg/colony, n=549, 52.3%), 0.4-0.6 (moderate, 11.2±5.8 kg/colony, n=249, 23.7%), and <0.4 (poor, 5.7±3.4 kg/colony, n=58, 5.5%). The substantial sample sizes across all categories, with the majority of observations(52.3%) falling within favorable conditions, demonstrate the practical utility and statistical reliability of these thresholds for beekeeping management decisions. AFI components exhibit nonlinear responses. Precipitation score (P_score) achieves maximum(1.0) under zero precipitation, declining to 0.6-0.8 for light rain (1-3 mm) and 0.1 for heavy rainfall (>5 mm), aligning with studies documenting rainfall’s mechanical flower damage, nectar dilution, and foraging suppression (Lee et al., 2020b). Diurnal temperature range (DT_score) peaks at moderate values (~12℃), reflecting requirements for circadian rhythm maintenance while avoiding thermal stress(Kim et al., 2020a). AFI advances beyond traditional Growing Degree Days (GDD) by incorporating multiple environmental factors. While recent North American Spring Index or European Unified Phenological Models include multiple variables, they employ linear combinations failing to capture nonlinear thresholds (Schwartz et al., 2006). AFI’s piecewise-defined scoring functions with empirically calibrated weights provide more flexible, biologically realistic frameworks. AFI addresses critical research gaps. The robustness of the AFI classification system is further supported by the adequate distribution of observations across all condition categories in the five-year dataset (2019-2025, N=1,050). Even in extreme conditions, including optimal scenarios (n=194) and severely constrained situations(n=58), there were sufficient observations to enable reliable statistical estimates. The substantial standard deviations observed across all AFI categories (SD ranging from 3.4 to 7.1 kg/colony) reflect the inherent variability in honey production, which is influenced by meteorological conditions, colony strength, apiary management practices, and local forage availability. Despite this variability, clear production trends emerged across AFI thresholds, demonstrating the dominant role of meteorological conditions as captured by the AFI in determining apicultural productivity. Earlier single-variable or simple combination studies achieved modest explanatory power (R²=0.40-0.55) (Kang and Lee, 2018). AFI’s systematic integration and nonlinear transformations substantially elevate predictive capacity (R²=0.76 training, R²=0.77 validation), approaching theoretical maxima given biological variability.

3. Stage-specific differentiation strategies: biological significance

Temporal stratification into three developmental stages-bud formation (-21 to -14 days), pre-flowering (-13 to -1 days), and flowering (0 to +12 days)-represents a conceptual advancement recognizing that meteorological impacts vary across developmental phases. Stage-specific approaches yielded substantially improved performance (R²=0.81) versus aggregated seasonal indices(R²=0.62).

The bud formation period constitutes a critical window when flower primordia undergo rapid cell division. AFI during this phase (AFI_bud) correlates strongly (r=0.68) with subsequent nectar volume and sugar concentration, suggesting that favorable early conditions establish physiological capacity that cannot be fully compensated later. This aligns with research demonstrating that carbohydrate allocation patterns during floral development exert lasting influences on reproductive investment (Tixier et al., 2013). The pre-flowering period shows diminished importance, reflected in lower regression coefficients for AFI_pre (β= -0.058, p=0.173). This counterintuitive negative association, though not significant, may reflect trade-offs where overly favorable conditions stimulate excessive vegetative growth competing with reproductive development, or indicate that moderate stress triggers adaptive nectar production responses. The flowering period emerges as most critical, with CI_flowering (bee activity index) dominating production models(β=0.929, p<0.001). This paramount importance reflects direct causal linkages between meteorological conditions, pollinator activity, and honey accumulation. Favorable flowering weather enables sustained foraging, efficient nectar collection, and rapid processing. CI’s correlation with production (r=0.872) substantially exceeds individual meteorological parameters or comprehensive AFI, highlighting pollinator behavior as the proximate mechanism linking environment to production.

Stage-specific modeling enables identification of vulnerable windows when adverse conditions disproportionately influence outcomes, providing essential tools for climate impact assessment and adaptive planning.

4. Practical applications and policy recommendations

The integrated STDD-AFI framework offers substantial utility across multiple stakeholder groups, from individual beekeeping operations to national agricultural planning and climate adaptation strategy.

For commercial beekeeping, the primary application lies in optimizing migratory strategies through accurate flowering forecasts. The framework’s predictive capacity, validated across 1,050 colony-year observations from five regions, provides statistically robust guidance for operational decisions. South Korea’s industry relies heavily on migration, with ~70% of operations moving colonies among multiple sites. STDD’s 3.8-day accuracy enables precise hive movement timing, ensuring colony arrival before flowering while avoiding premature deployment costs. Regional AFI forecasts provide quantitative route planning guidance, helping prioritize destinations with superior conditions while avoiding poor production areas. Regression-based production forecasting enables data-driven harvest expectations and marketing strategies. Early-season meteorological observations during bud formation and pre-flowering inform preliminary estimates months before harvest, facilitating buyer negotiations and financial planning. Progressive incorporation of flowering-period data allows iterative forecast refinement with increasing precision.

For research institutions and extension services, the framework provides foundations for accessible decision support tools. Web-based or mobile platforms integrating STDD-AFI models with real-time meteorological streams could deliver personalized recommendations based on geographic locations. Alert functions notifying flowering onset, quality ratings for upcoming periods, and comparative regional assessments would enhance operational efficiency. From policy perspectives, the framework addresses critical climate adaptation planning needs. Government agencies can utilize projections to inform targeted support programs, directing resources toward regions with elevated climate-related risks. Subsidized insurance programs could employ AFI-based triggers for eligibility and payouts, providing financial protection against meteorological production failures while incentivizing climate-adaptive practices.

Long-term monitoring programs integrating standardized phenological observations with systematic meteorological collection provide essential foundations for adaptive management. The National Institute of Agricultural Sciences’ coordinated survey (77 fixed locations, 144 migratory operations with GPS since 2019) exemplifies necessary data infrastructure. Sustained commitment with periodic model recalibration ensures decision support relevance as conditions evolve. Realizing full potential requires complementary capacity building and technology transfer investments. Many small-scale traditional beekeepers lack internet connectivity or technical literacy for sophisticated tools. Extension programs must employ diverse strategies-printed materials, radio broadcasts, workshops, peer learning-ensuring equitable forecasting access. Demonstration projects showcasing model-guided successes can build trust and encourage broader data-driven adoption.

5. Research limitations and future directions

While representing substantial advancement in integrated meteorological-phenological-production modeling, several limitations warrant acknowledgment and suggest future priorities.

First, temporal scope (2006-2017) may inadequately capture recent acceleration in climate change impacts and extreme event frequency. The past decade witnessed unprecedented warming trends, increased late-spring frosts, and altered precipitation patterns potentially shifting statistical relationships underpinning STDD. Periodic recalibration incorporating recent observations, particularly extreme years, would enhance robustness. Ideally, extending calibration datasets to 20-25 years would adequately represent interannual climate variability including rare high-impact events.

Second, spatial resolution (500 m × 500 m) constrains applicability to individual apiary-scale decisions. South Korea’s mountainous terrain generates pronounced microclimatic heterogeneity inadequately resolved at this scale. Downscaling approaches incorporating high-resolution elevation models and land cover could refine predictions, though validation requires expanded monitoring networks including distributed automated sensors.

Third, the framework focuses exclusively on meteorological drivers while neglecting non-meteorological influences. Despite achieving strong predictive performance (R²=0.76) with adequate sample sizes across all regions (n=210 per region), apiary management practices (queen quality, hive strength, disease, supplemental feeding, varroa mite infestation) substantially impact productivity independent of environment and likely contribute to the observed standard deviations in production outcomes (SD=3.4-10.1 kg). Landscape factors (alternative forage, pesticide exposure, agricultural proximity) may modulate meteorological responses. Integrating such variables requires interdisciplinary data collection encompassing management surveys, landscape characterization, and veterinary assessments.

Fourth, climate adaptation requires projecting future conditions under evolving regimes. Empirically calibrated models assume stationarity in physiological responses and phenological cue sensitivities, which may not hold as organisms adapt through plasticity or evolution. Mechanistic process-based models grounded in biophysical principles offer greater confidence for extrapolation beyond observed ranges. Hybrid approaches combining mechanistic core processes with empirical calibration, potentially employing machine learning (random forests, gradient boosting, neural networks), could capture complex nonlinear relationships while maintaining interpretability.

Fifth, the framework considers R. pseudoacacia in isolation, neglecting interactions with other spring-flowering species competing for pollinators or providing supplementary forage. In mixed forests or agricultural landscapes with diverse phenologies, production depends on synchronization with competitors, sequential forage availability, and landscape-scale resource abundance. Expanding frameworks to incorporate community-level phenological dynamics would enhance applicability.

Future research should focus on several key directions. First, incorporating climate scenario projections (SSP1-2.6, SSP2-4.5, and SSP5-8.5) will enable assessment of long-term viability, accounting for both gradual climatic trends and interannual variability, as well as the frequency of extreme events. Second, developing a mechanistic understanding of key physiological processes through controlled experiments is essential, particularly those examining thermal response curves, precipitation effects on nectar sugar concentration, and wind thresholds that influence foraging cessation. Third, applying machine learning approaches offers potential for integrating diverse datasets and producing probabilistic forecasts with improved precision. Fourth, expanding the model scope to include additional honey flow periods, such as summer chestnut and autumn aster seasons, would allow for more comprehensive annual yield forecasting. Finally, future efforts should investigate how climate factors affect honey quality, including variations in floral composition, color, antioxidant levels, and antimicrobial properties, beyond mere production quantity.