The 2019–2020 Australian bushfires released approximately 715 million tonnes of carbon dioxide (CO₂), equivalent to over two-thirds of Australia’s annual fossil fuel emissions (van der Velde et al., 2021). This substantial release of greenhouse gases intensified global warming and raised atmospheric CO₂ levels. The wildfires also resulted in the loss of 21% of Australian temperate forests, significantly reducing the carbon storage capacity of these ecosystems (Shiraishi and Hirata, 2021). Over the past decades, southeastern Australia has experienced a 30% increase in fire weather severity since 1950, driven by record-breaking temperatures and prolonged droughts. Future climate scenarios projected a worsening trend, with the number of extreme fire weather days expected to double by the end of the century (Abram et al., 2021). These developments underscored the urgent need to monitor vegetation recovery to address long-term ecological impacts (Bradstock, 2010).

Previous studies that monitored vegetation recovery after the wildfire relied on satellite datasets with low spatial resolution or low temporal resolution. While these datasets effectively captured large-scale and long-term changes, they faced limitations in identifying fine-scale recovery patterns. For instance, Peña-Molina et al. (2024) analyzed changes in the Normalized Difference Vegetation Index (NDVI) and Normalized Burn Ratio (NBR) for fire-affected regions in Australia using Landsat data. Although their analysis successfully identified broad recovery patterns, it highlighted significant limitations in detecting delayed or failed recovery areas during early regrowth stages. These limitations stemmed from the low temporal coverage of Landsat, which was insufficient to capture short-term dynamics.

This study aimed to address the spatial and temporal resolution gaps inherent in traditional satellite products by generating highresolution datasets through image fusion. By integrating Moderate Resolution Imaging Spectroradiometer (MODIS) MCD43A4 products with Landsat 8/9 Nadir Bidirectional Reflectance Distribution (NBAR) data, we created daily gap-free imagery at a 30-meter resolution, overcoming the limitations of individual datasets. To evaluate the performance of the fused dataset, we assessed its ability to detect vegetation recovery patterns in Tumbarumba, a wildfire-affected region in southeastern Australia. Additionally, we compared vegetation recovery insights derived from the fusion dataset with carbon flux measurements obtained from a 70-meter eddy covariance tower in Tumbarumba (Leuning et al., 2005). This study specifically addressed two research questions: whether fusion imagery outperforms traditional satellite products in monitoring wildfire recovery and how the temporal recovery pattern of vegetation in the target site appears.

2.1. Study Area and Period

Tumbarumba, located within the Bago State Forest in southeastern New South Wales, Australia, was selected as the study area (centered at Global Positioning System [GPS] coordinates with latitude: –35.6566° and longitude: 148.1517). The research site covered a 10 × 10 km area around these coordinates (Fig. 1). The forest was classified as a wet sclerophyll forest and was dominated by Alpine Ash (Eucalyptus delegatensis), with an average tree height of 40 meters. The site was situated at an altitude of 1,200 meters and receives an annual mean precipitation of approximately 1,000 mm.

Figure 1. Burned area map from Google Earth Engine provided by the New South Wales (NSW) government (NSW Department of Planning and Environment, 2020). The fire intensity was classified into six levels by color, with the study site marked by a black rectangle.

Figure 2. The entire process of image fusion. MODIS was used as coarse images, and BRDF correction and gap filling were performed using Landsat as fine images.

The study period was 10 years, from 2014 to 2023, covering both before and after the wildfire conditions. On December 31, 2019, a moderate-severity bushfire burned through the Bago State Forest, consuming the understory vegetation while leaving the canopy intact. Alpine Ash experienced nearly 100% mortality, while Mountain Gum (Eucalyptus dalrympleana) resprouted with epicormic growth(Keith et al., 2009). Regrowth after the wildfire included Alpine Ash seedlings and resprouting of other eucalyptus species such as Mountain Gum and Peppermint (Ozflux, 2022). A 10-month gap in data collection occurred immediately following the wildfire.

2.2. Data Availability

2.2.1. Satelite Data

The MODIS provided daily, high-temporal resolution satellite data through NASA’s Earth Observation Program (NASA, 2023). This study utilized the MCD43A4 product, which applied the Bidirectional Reflectance Distribution Function (BRDF) correction to normalize surface reflectance values for variations in solar illumination and sensor viewing geometry. The RossThick- LiSparse-Reciprocal model achieved BRDF correction in the MCD43A4 product. This semi-empirical model combined isotropic, volumetric, and geometric scattering kernels to account for angular effects in surface reflectance (Schaaf et al., 2002). The correction adjusted surface reflectance values to nadir-viewing conditions, generating the NBAR to remove directional biases inherent in the observation geometry. The MCD43A4 dataset aggregates multi-angle observations over 16 days to construct BRDF models for each pixel, allowing for the computation of NBAR under consistent viewing conditions. This product, available from the NASA Land Processes Distributed Active Archive Center (LP DAAC), delivers reliable daily surface reflectance data at a spatial resolution of 500 meters. All available MODIS NBAR datasets over the study area during the 10 years were utilized in the fusion process. To ensure the temporal consistency of MODIS NBAR, we used both high-quality BRDF inversion and magnitude BRDF inversion products.

This study utilized Landsat 8 and Landsat 9 imagery, which provided high spatial resolution satellite data through the Earth observation program managed by the United States Geological Survey (USGS) (USGS, 2023). The dataset includes Level 1 and Level 2 products, each offering specific features tailored to different analytical needs.

Level 1 data were radiometrically and geometrically corrected, which ensures pixel alignment with global geographic reference systems and enables precise spatial analyses. Additionally, Level 1 products contained observation geometry information, such as solar and satellite view angles, which are crucial for BRDF corrections or angular normalization of reflectance values (Roy et al., 2014). This study utilized the angle information from Level 1 data to account for geometric distortions and angular dependencies when analyzing surface reflectance.

Level 2 data provided surface reflectance values adjusted for the effects of aerosols, water vapor, and ozone. This processing ensured that reflectance measurements represent true surface properties, which makes the data ideal for calculating vegetation indices like NDVI and monitoring temporal changes in vegetation recovery (USGS, 2023). This study utilized angle information from Level 1 data and surface reflectance values from Level 2 data across the visible (blue, green, red) and near-infrared (NIR) bands to generate NBAR products. These NBAR datasets were used to analyze vegetation recovery following the wildfire.

Over the 10-year study period, a total of 115 Landsat images were used, averaging approximately 11.5 images per year and 0.96 images per month. Due to Landsat’s 16-day revisit cycle, the availability of images varied depending on cloud cover and acquisition schedules, leading to months with no usable images.

2.2.2. Image Fusion

To address the challenges of limited temporal and spatial resolution in satellite imagery, a workflow integrating image fusion and gapfilling methods was employed. The purpose of this process was to generate daily, gap-free imagery with a 30-meter spatial resolution, combining the high temporal resolution of MODIS and the fine spatial resolution of Landsat to facilitate detailed environmental analysis. The workflow comprised three main stages: preprocessing, fusion, and gap filling.

The preprocessing stage was designed to prepare MODIS and Landsat data for use in the image fusion process. For MODIS, the Neighborhood Similar Pixel Interpolator (NSPI) algorithm was applied to reconstruct missing data by identifying spectrally similar neighboring pixels, resulting in daily, gap-free coarseresolution images (Zhu et al., 2012). Landsat data underwent preprocessing steps that included BRDF correction to normalize surface reflectance values. This process utilized angle data from Landsat Level 1 products and surface reflectance data from Level 2 products. BRDF correction adjusted Landsat reflectance values to nadir-viewing conditions using MODIS BRDF model parameters, thereby aligning the datasets for the fusion process. This adjustment normalized reflectance to nadir-viewing and solar-normalized conditions by accounting for angular effects caused by varying observation and solar geometries (Roy et al., 2008; 2016). Following BRDF correction, Landsat data were cropped to the study area, and cloud-contaminated pixels were detected and masked using the Fmask algorithm. Fmask generated masking files by identifying cloud, and cloud shadow through spectral-contextual methods (Qiu et al., 2019).

The fusion stage synthesized the gap-filled MODIS and preprocessed Landsat data, focusing specifically on the NIR and visible bands (red, blue, green), using an enhanced FSDAF that incorporates sub-pixel class fraction change information (SFSDAF) algorithm. SFSDAF predicts high-resolution imagery at 30 meters by modeling temporal reflectance variations from MODIS and spatial heterogeneity from Landsat. This algorithm incorporated sub-pixel class fraction changes, capturing both gradual phenological shifts and abrupt land cover changes (Li et al., 2020). To ensure quality in the fusion process, only Landsat images with less than 50% cloud cover were included in the dataset, minimizing contamination from cloud-obscured regions. However, the influence of clouds was not eliminated.

In the final gapfilling stage applied the NSPI algorithm to address residual gaps in the fused dataset. NSPI interpolates missing values based on spectral similarity and spatial proximity while utilizing the masking files generated during the Fmask process to refine gap detection and filling. This step specifically targeted areas with cloud-induced missing data, ensuring the production of a 30-meter daily gap-free dataset.

2.2.3. Eddy Covariance Data

To compare satellite data with flux tower observations for analyzing vegetation recovery, this study utilized Tumbarumba flux tower data processed to Level 6 (L6) format. The original Level 3 (L3) data, provided by the OzFlux network, were further refined using the REddyProcWeb online tool, developed by the Max Planck Institute for Biogeochemistry (Wutzler et al, 2018). This tool applied widely recognized methodologies for partitioning Net Ecosystem Exchange (NEE) into Gross Primary Production (GPP) and ecosystem respiration, following the approach described by Reichstein et al. (2005). The partitioning process leverages nighttime NEE data and temperature sensitivity to estimate ecosystem respiration, from which GPP is calculated as the residual flux. The REddyProcWeb tool also implemented gap-filling techniques to address missing data in the flux tower observations.

2.3. NDVI

To analyze vegetation recovery and dynamics, this study utilized the NDVI. NDVI is a widely used remote sensing index for assessing vegetation health and greenness, calculated using the formula:

NDVI is a widely used index that quantifies vegetation health by comparing the strong reflection of NIR light with the absorption of red light by chlorophyll for photosynthesis. NDVI values typically range from 0 to 1, where higher values indicate denser and healthier vegetation, while lower values correspond to sparse vegetation, bare soil, or water (Tucker, 1979).

3.1. Comparing Spatio-Temporal Resolution

The fusion algorithm generated daily gap-free images at 30-meter spatial resolution. Compared to MODIS, which provided daily data at a coarse resolution of 500 meters, the fusion dataset provided spatial detail, allowing finer landscape features to be identified. In contrast, Landsat imagery, while offering high spatial resolution, was limited by its 15-day revisit interval. Landsat did not provide data for dates between revisits, and the availability of usable images further decreased when cloud cover occurred. In contrast, the fusion dataset provided consistent daily gap-free images across all periods.

3.2. Monitoring Vegetation Recovery

3.2.1. Temporal Dynamics of NDVI

Monthly median NDVI values for the Tumbarumba were analyzed at 480-meter and 30-meter resolutions using fusion images, MODIS, and Landsat data (Fig. 4). At 480 meters, both the fusion dataset and MODIS data displayed a similar temporal pattern, including a sharp NDVI decline following the wildfire on December 31, 2019, and subsequent recovery. Two years after the wildfire, the fusion images recorded NDVI recovery to 98% of the levels observed before the wildfire, while MODIS indicated a 95% recovery (Fig. 4a). At 30 meters, the fusion dataset revealed that two years after the wildfire, the average NDVI had recovered to 83% of the levels prior to the wildfire (Fig. 4b). Additionally, the fusion images demonstrated better consistency in capturing seasonal patterns compared to Landsat.

Figure 4. NDVI time series at the Tumbarumba tower site, comparing the fusion image with (a) MODIS at 480 m resolution and (b) Landsat at 30 m resolution.

Spatial patterns of recovery were analyzed using 10 × 10 km NDVI maps generated for each January over 10 years (Fig. 5). The January 2020 map indicated a marked decline in NDVI following the December 31, 2019 wildfire, reflecting significant vegetation loss. In the subsequent years, NDVI values showed a consistent increasing trend, though the degree of increase varied across different locations.

Figure 5. 10 x 10 km NDVI maps of Tumbarumba for each January over 10 years, derived from (a) MODIS images and (b) fusion images.

3.2.2. Temporal Dynamics of Eddy Covariance Observations

The wildfire caused changes in GPP and NEE values before and after the event. GPP decreased from an average of 8.19 μmol CO₂ m^–2 s^–1 during 2014–2019 to 4.61 μmol CO₂m^–2 s^–1 during 2020– 2023, representing a 43.74% decline (Fig. 6a). In contrast, the NEE shifted from a mean of –2.76 μmol CO₂ m^–2 s^–1 before the wildfire to –1.19 μmol CO₂m^–2 s^–1 after the wildfire, representing a 56.85% increase (Fig. 6b).

Figure 6. Comparison of before and after wildfires. (a) GPP mean and (b) NEE mean.

4.1. Patterns of Vegetation Recovery After Wildfire

NDVI calculated from the fusion dataset demonstrated a rapid increase in vegetation greenness in the Tumbarumba after the December 31, 2019 wildfire. Within two years, NDVI values recovered to 98% of the values before the wildfire at 480-meter resolution and 83% at 30-meter resolution. This recovery showed a significant rebound in vegetation cover. These findings aligned with previous studies, which showed that ecosystems affected by wildfires often experience rapid regrowth due to adaptive plant strategies and favorable conditions following the wildfire.

However, flux tower data revealed contrasting trends. GPP declined by 43.74% from 8.19 μmol CO₂ m^–2 s^–1 during 2014–2019 to 4.61 μmol CO₂ m^–2 s^–1 during 2020–2023, indicating reduced photosynthetic capacity. Similarly, NEE shifted from –2.76 μmol CO₂m^–2 s^–1 before the fire to –1.19 μmol CO₂m^–2 s^–1 after, reflecting diminished carbon sequestration capacity (Fig. 6). These results highlighted that while NDVI effectively captures vegetation greenness, it does not represent physiological recovery processes such as photosynthesis or carbon exchange. Observing both NDVI with flux tower data provided a more nuanced understanding of vegetation recovery after the fire.

4.2. Importance of High Spatio-Temporal Resolution in Vegetation Monitoring

High spatial resolution played a critical role in detecting detailed variations in vegetation recovery dynamics that were not observable at coarser resolutions. At 480-meter resolution, recovery trajectories appeared uniform, potentially masking smaller-scale variations in vegetation response to the wildfire. This limitation could result in overlooking localized areas of delayed recovery or persistent vegetation loss, which are crucial for understanding vegetation resilience. In contrast, the 30-meter resolution data revealed spatial heterogeneity in NDVI recovery, identifying areas with slower regrowth and sustained vegetation loss (Fig. 5). These findings suggest that higher spatial resolution provides key insights into site-specific variations, which can inform targeted restoration strategies.

Temporal resolution was equally critical for capturing dynamic changes in vegetation over time. The fusion dataset, offering daily observations, captured consistent seasonal patterns and provided detailed insights into short-term vegetation dynamics. Landsat data, constrained by its 15-day revisit interval, often missed these nuances, particularly during early recovery phases. During these periods, rapid changes in vegetation indices occurred due to regrowth, seasonal variations, and environmental factors. The daily temporal resolution allowed the detection of recovery spurts during optimal growing seasons, which are crucial for understanding recovery trajectories and assessing intervention effectiveness (Fig. 3).

Figure 3. Resolution and temporal availability of MODIS, Landsat, and fusion images from May 20, 2022, to August 16, 2022.

Despite these advantages, maintaining the quality of high temporal resolution datasets posed challenges. For example, as observed in Fig. 4, outliers were detected in May 2014 and June 2023, likely caused by cloud interference or the preprocessing process. Cloud contamination often caused missing data outliers that were not fully removed, disrupting observed trends. Detecting these outliers remained a difficult task, particularly in dynamic environments after the wildfire where rapid changes may mimic anomalies. While the dataset effectively captured seasonal patterns, the lack of applied temporal smoothing occasionally led to inconsistencies in data continuity. Implementing temporal smoothing methods in future studies can help address these inconsistencies, ensuring that seasonal variations are accurately represented without oversimplifying the complex recovery dynamics. Finally, in terms of data availability, this study has potential biases. Landsat data tend to have a large number of data available mainly in winter and as little as half in summer (Fig. 7). This seems to be due to an increase in the amount of cloud in the summer, resulting in an increase in the time interval. These changes in temporal resolution can lead to difficulties in detecting rapid changes in the period (Wang and Atkinson, 2018). Further use of satellite imagery with lower time inverter in the future will minimize this problem. These improvements are essential for maximizing the reliability and interpretability of high-resolution fusion datasets in long-term vegetation monitoring.

Figure 7. Number of available Landsat data per month during the 10 years (2014–2023).

This study emphasized the effectiveness of high-resolution image fusion in monitoring vegetation recovery after wildfires. The fused dataset successfully captured fine-scale spatial variability and temporal dynamics that could not be fully addressed by MODIS and Landsat products. NDVI observations revealed significant differences in vegetation recovery between resolutions. At 480-meter resolution, recovery appeared nearly complete, with NDVI values reaching 98% of the levels observed prior to the wildfire within two years. In contrast, the 30-meter resolution data revealed a lower recovery of 83%, highlighting localized variations in recovery rates that were otherwise masked at coarser resolutions. This demonstrates that high-resolution satellite imagery enables the identification of small-scale recovery units and localized patterns that are crucial for understanding heterogeneous ecosystem responses to wildfires. Ground-based measurements, such as flux tower data, provided additional insights into the recovery of photosynthesis and carbon sequestration. Despite these advancements, challenges such as outliers and temporal inconsistencies continue to persist. Improving gapfilling and temporal smoothing algorithms will be critical to enhancing the reliability of fusion datasets. This study underscores not only the importance of high-resolution image fusion for vegetation recovery monitoring but also the critical role of integrating diverse datasets in enhancing our understanding of vegetation recovery dynamics.

This research was supported by the Technology Development Project for Creation and Management of Ecosystem-based Carbon Sinks (202300218237) through the Korea Environmental Industry & Technology Institute (KEITI) funded by the Ministry of Environment (MOE).

No potential conflict of interest relevant to this article was reported.