The YS, located between China and Korea, is a semi-enclosed marginal sea in the Pacific Northwest with an average depth of 44 m (Jin et al., 2013). The YS is predominantly affected by aerosols, as it is near the northeastern part of China where there is significant industrial activity. As a result, this area is influenced by atmospheric nitrogen deposition and it exhibits high primary productivity (Yoon et al., 2021). In addition, the YS is influenced by river inflows containing suspended sediment. This study examined the influence of aerosols on phytoplankton growth in an area that minimizes the effects of river discharge and suspended sediment. The central region of the YS (Central Yellow Sea [CYS]; 121.5–125.5°E, 34–37°N) was selected (Figure 1) (Shi et al., 2017; Fu et al., 2018) because it was the region where spring blooms mainly occurred in the YS (Xu et al., 2013). In addition, the sources of land-based AOD in the study area did not make a substantial contribution to atmospheric aerosols in the CYS, as per the ocean currents and wind directions. Hence help intensively investigate the effect of aerosols on the CHL. According to Xu et al. (2013), spring blooms in the CYS were defined as a CHL higher than 1.2 mg m^-3.

According to previous studies, among the ML techniques, the RF model simulates fairly well the complex effect of various factors on CHL (Liu et al., 2013; Liu and Wu, 2017; Park et al., 2019). Thus, we performed a regression analysis to evaluate the CHL response to AOD and other physical forcing variables (SST, MLD, WD, and TP) of the RF model. The RF model, an ensemble learning method, can perform classification and regression analyses to evaluate the CHL response to AOD and other physical forcing variables (SST, MLD, WD, and TP). Decision tree-based algorithms, such as RF, can measure the variable importance (VI) of the features used in model development (Breiman, 2001; Ali et al., 2012). After training, this score was calculated for each feature and the result was normalized such that the total sum of importance was 1. Thus, the VI, which indicates the effect of the predictors used in the RF model on CHL, could be evaluated (Cutler et al., 2007; Park et al., 2020; Pham et al., 2021).

Generally, the RF model sets the number of single decision trees (N_trees) and the number of features (M_try) to generate the best performance in each tree node. For regression, M_try was set to the maximum number of features and for classification, it was set to the square root of the number of features (Geurts et al., 2006). Then, for the feature selected at each node, the node was divided into two child nodes based on the optimal partitioning condition (i.e., a small mean squared error, MSE) and a leaf node was calculated (Park et al., 2020). Finally, the results produced by each tree were ensembled on average to draw a conclusion. The critical parameters considered in the RF model were N_trees and M_try. N_trees was set to 50 through an empirical trial and M_try was set to 6 according to the convention.

Six predictors (variables) were used as the inputs for the RF model. The predictors used were the AOD, physical forcing variables (SST, MLD, WD, and TP), and climatology of the CHL (CHL_c) during spring (March 1 to March 31) from 2003 to 2021. As inputs for the RF model, we used CHL_c because it is considered the factor illustrating the environmental variables (Park et al., 2020). It is worth noting that PAR was not used; instead, we used AOD as input because both have directly inverse relations (Li et al., 2020). The target was the MODIS-Aqua CHL (CHL_MODIS). The significance test of the predictors was conducted using the Boruta algorithm and all the data were normalized. In addition, for CHL_c and CHL_MODIS, which have skewed data, a logarithmic transformation was performed to approximate a normal distribution (Gregor et al., 2017; Hahm et al., 2019).

In Section 3.2, the performance of the RF models was evaluated using various statistical metrics, such as the determined coefficient (R²), root mean square error (RMSE), mean absolute error (MAE), and MSE (Popoola et al., 2019; Park et al., 2020). The R² indicates how well the estimated linear model fits the given data in statistics. RMSE is a measure commonly used when dealing with the difference between the in-situ data (X) and predicted value (Y). MAE means errors between paired observations expressing the same phenomenon. In other words, MAE is the mean absolute difference between the X and Y. MSE shows the mean of the squares of the error. The closer the MSE is to 0, the higher the accuracy is because the guessed value is closer to the original (Chicco et al., 2021). In addition, using the RF model, we determined whether CHL could increase when only AOD was additionally supplied to environmental conditions in March 2020 (Section 3.3). Subsequently, the quantitative contribution of each factor, including AOD to CHL in the RF model, was evaluated using a partial dependence (PD) plot (Section 3.4). A partial dependence (PD) plot depicts the relationship between input variables and predictions (Friedman. 2001). These results show how the predictions partially depend on the values of the input variables of interest. This method can be used to interpret a trained model that can identify linear, monotonic, or more complex relationships between predictors and targets (Molnar, 2020). Furthermore, the PD plot quantitatively describes the response probabilities of the predictors (Friedman, 2001; Greenwell, 2017; Zeng et al., 2017). In other words, the PD plot describes the dependence between the target and predictors of interest, marginalizing over the value of all other predictors. We considered that this PD analysis could represent the quantitative effect on the absolute contribution of each physical variable to CHL, that is, the amount of change in CHL (Park et al., 2020). PD analysis was performed as follows:

where f^(x)x={x1x2⋯xp} represents the number of features in the RF model (p = 6 in this study). In addition, when we divide x by the set of selected X_s and complement (X_c ) the response variable in the PD plot for X_a is defined by Eq. (1). In other words, X_s is the feature (i.e., variables and inputs) to be calculated through the model and X_Cis all variables except X_s used feature. p_c(X_c) represents the marginal probability density of X_c.X_s is set to one variable, and X_c is the average value that the remaining inputs are affected by the change of X_s . This means that the effects of all the other predictors in the model were average. For more details, refer to Park et al. (2020). Eq. (1) can be evaluated using the training set, as follows:

where X_i,c (i=1, 2, ..., n) represents the value of X_c in the training samples.

The methodology of this study was based on that of Park et al. (2020) (Park et al., 2019). However, unlike the referenced papers on the RF approach, in this study, regression analysis was used among the RF model features to understand the effect of each factor on CHL. In addition, after quantitatively evaluating the impact of each element on spring blooms, the effect was formulated as an equation. Accordingly, this study attempted to extract critical information on how AOD was quantitatively related to spring blooms during the COVID-19 lockdown period.

First, we determined how much the aerosol decreased during the COVID-19 lockdown period compared to the other years; therefore, we compared AOD and CHL in March from 2019 to 2021 (Figure 1). During the COVID-19 lockdown period in 2020, local patches were inferred to be small spring blooms (Figure 1E). Especially, it was confirmed that the levels of both data in March 2020, during the COVID-19 lockdown period, decreased significantly compared to those of the other years (Figures 1B, E). In 2019 and 2021, we observed that regions with high AOD distribution showed high CHL levels, confirming that both variables had similar spatial distributions (Figures 1A, C, D, F). From Figure 1, it can be inferred that the abnormally low CHL levels in 2020 was affected by low AOD. We also investigated the changes in other factors to determine whether the effect of reduced aerosol weakened spring blooms during the COVID-19 lockdown period.

Figure 2 shows that AOD and CHL in March 2020 were lower than the average values from March 2003 to 2021, by approximately 50% (from 0.50 to 0.26) for AOD and approximately 15% (from 2.18 mg m^-3 to 1.90 mg m^-3) for CHL (Figures 2A, F). In contrast, Figure 2B shows that SST increased by approximately 20% from 7.61°C to 9.34°C during the same period. This warmer SST is likely due to the reduced amounts of aerosols, which might have caused the increase in SST due to higher incoming solar radiation (Yue et al., 2011; Kripalani et al., 2022). Accordingly, an increased SST can induce strengthened thermal stratification and a shallower MLD (Tan and Shi, 2012a; Lu et al., 2021). Consequently, the MLD in 2020 became 25% more shallow, from 44.37 m to 33.09 m (Figure 2C). In addition, phytoplankton biomass generally decreases as the MLD becomes more shallow (Lu et al., 2021). This was also seen for the CHL in 2020 and 2021. Additionally, the WD indicated that the airflow from the atmosphere to the ocean was slightly weakened by 13% (from -2.13 × 10^-6 sec^-1 to -2.39 × 10^-6 sec^-1), but this value was not significantly different than that of other years. TP also decreased by 45%, from 7.38 × 10^-5 mm to 3.92 × 10^-5 mm, which is presumed to be an effect of aerosol reduction in humide environmental conditions (Figures 2D, E) (Khain et al., 2008).

In the case of 2021, many factors showed a significant change compared with the average values from March 2003 to 2021. AOD showed a sharp increase compared with the same period, and SST showed a similar increasing pattern in 2020. Therefore, the thermal stratification was strengthened owing to the increased SST, and the MLD became shallow, as in the COVID-19 lockdown period. In contrast, the WD and TP in 2021 showed much higher levels than the mean values throughout the period, and the increased WD strengthened the flow from the atmosphere into the ocean. The increase in TP may have been because of increased aerosols, and atmospheric aerosols may have possibly directly entered the ocean through precipitation. Consequently, in 2021, with an increase in aerosols, more aerosols could flow into the sea from the atmosphere due to the physical environment’s conditions. We confirmed that 2021 showed a higher CHL than 2020 under different physical environmental conditions during the COVID-19 lockdown period (e.g., increased SST and decreased MLD). Accordingly, the spring blooms in the study area were mainly controlled by AOD changes, as shown in Figures 1, 2. Ultimately, we focused on the qualitative and quantitative effects of aerosol changes that significantly influenced spring blooms.

To identify the effect of each factor on CHL, we developed a RF model using the above-mentioned data. The performance of the training set with 80% of the total data (n = 94,480) showed R² = 0.98, RMSE = 0.22 mg m^-3, MAE = 0.07 mg m^-3, and MSE = 0.05 mg m^-3 with high accuracy and low RMSE (Figure 3A; Table 2). It is worth noting that the splitting ratio of 80% to 20% for training and validation dataset shows slightly better results compared to other splitting ratios. According to Despotovic et al. (2016), relative RMSE (RRMSE) evaluates model construction accuracy and is considered good when 10 % ≤ RRMSE < 20 % (i.e., 10.49 % indicates a good level), excellent if RRMSE < 10 %, fair if 20 % ≤ RRMSE < 30 %, and poor if 30 % ≤ RRMSE (Jamieson et al., 1991; Heinemann et al., 2012; Li et al., 2013). Hence, the training set almost coincided with CHL_MODIS and the reconstructed CHL derived from the RF model (CHL_RF) (Figure 3A). The performance of the validation set with 20 % of the total data (n = 23,620) had accuracy with R^{2 =} 0.85, RMSE ⁼ 0.54 mg m^-3, MAE = 0.17 mg m^-3, and MSE = 0.29 mg m^-3, and RRMSE = 25.37 %. Although the validation set had a lower R² and higher RMSE than those of the training set, thus, the reconstructed result was significant (Figure 3B; Table 2). Additionally, we determined the relative contribution of each predictor to the RF model through the VI (Figure 3C). The relative contributions of CHL_c, TP, SST, WD, MLD, and AOD were 46.10 %, 14.42 %, 12.87 %, 9.64 %, 8.71 %, and 8.27 %, respectively. In addition, the standard deviations of CHL_c, TP, SST, WD, MLD, and AOD were ±0.29 %, ±0.47 %, ±0.44 %, ±0.44 %, ±0.49 %, and ±0.39 %, respectively. CHLc, the average CHL from March 2003 to 2021, had the most significant effect on the CHL prediction because it reflects the CHL trend. Hence, CHL_c was excluded from the physical forcing variables that directly influenced CHL variation.

To determine whether the RF model could simulate the regional pattern, we compared the spatial distributions of CHL_MODIS and CHL_RF (Figure 4). The locations and intensities of the spring blooms were also closely simulated each year. The CHL_MODIS was 2.09 mg m^-3, 1.89 mg m^-3, and 2.71 mg m^-3, and CHL_RF was 2.01 mg m^-3, 1.82 mg m^-3, and 2.66 mg m^-3 in 2019, 2020, and 2021, respectively. If we estimated CHL_RF without CHL_c, the CHL_RF predicted an overestimate (not shown). This result may be because CHL_c reflects the CHL information trend. In particular, spring blooms in 2020, which occurred on a small scale and were weak compared to those in other years, were simulated reasonably well. Therefore, it was confirmed that the RF model was well-reconstructed based on the intensity and spatial pattern of CHL (Figures 4B, E).

The established RF model simulated the spatial pattern and CHL value using six factors. The five factors that directly affect CHL had slight differences, but their relative contributions to the RF model were almost similar. Previous studies have shown that CHL is positively related to SST (Zhao et al., 2019) and negatively related to the MLD and TP (Liu et al., 2013; Lu et al., 2021). The influence of each factor (SST, MLD, WD, and TP), which changed as AOD decreased, was discussed in Section 3.1. Therefore, we concentrated more on the spring blooms based on CHL affected by AOD, which changed the most during the COVID-19 lockdown.

We focused on the relative effect of AOD on spring blooms, which contributed the least to the RF model but showed the most significant change in 2020 (Section 3.2). This section describes the sensitivity experiments for CHL_RF to AOD variations under artificially changing aerosol levels (Figures 4, 5). CHL_RF satisfactorily represents the spring blooms that occur each year. However, spring blooms had lessened in the COVID-19 lockdown period of 2020 compared to that in other years (Figure 4E). Therefore, we attempted to confirm the change in spring blooms by artificially increasing and decreasing aerosols during the COVID-19 lockdown period through a sensitivity experiment using the RF model.

Two sensitivity experiment conditions were constructed to investigate changes in CHL variation by adjusting AOD. The two conditions are with AOD ± two sigma (2 σ ± 95%), according to the empirical rule (three sigma rule) of statistics; where 2 σ is derived by assuming the Gaussian fitting for the distribution of AOD climatology. Furthermore, the sensitivity experiments are conducted with AOD + σ and AOD − σ while all other predictors remain unchanged. In other words, when only the AOD either increased or decreased, all other predictors except AOD are constant. The CHL changes to high (CHL_HAOD) and low (CHL_LAOD) AOD were estimated using the following equations:

The CHL_HAOD increased by approximately 95% in AOD and showed an increased CHL in most areas when compared with the existing CHL_RF (Figures 5A–C). When aerosols were artificially increased during the COVID-19 lockdown period (in March 2020), phytoplankton biomass increased rapidly in most regions compared with the biomass values from other years (Figure 5B). In particular, although differences occur depending on the region, CHL_RF increased by about 60% or more on average in 2019, 2020, and 2021. Based on this result, if the AOD had additional supplies, it could have generated intense spring blooms in March 2020. The CHL_LAOD decreased by approximately 95% in AOD and also exhibited a decreased CHL in most areas when compared with CHL_RF (Figures 5D–F). Furthermore, Figures 5D, 5F show that the region had significantly lower CHL when the AOD was decreased by approximately 95%. However, the CHL_LAOD in 2020 did not represent a significant decrease in most areas when compared to CHL_RF because the AOD in March 2020 was already reduced by approximately 50% (0.50 to 0.26) (Figure 5E). This result indicates that the effect of AOD on CHL during the COVID-19 lockdown period was analogous to that of CHL_LAOD. This also proves that the experimental results suitably simulated the CHL response to the actual change in the aerosols in the environment.

Therefore, in Section 3.3, we presented the experimental results for predicting CHL changes according to real environmental conditions through CHL reactions, based on artificial AOD fluctuations. These results suggest that AOD (i.e., atmospheric aerosol deposition) significantly affects spring blooms in the CYS.

In Sections 3.2 and 3.3, we calculated the relative contribution of each factor to spring blooms. Relative contribution refers to the relative value that each factor can influence when various physical factors affect the CHL in combination. However, in Section 3.4, the quantitative contribution was considered as the absolute impact of each variable on CHL individually. This means that the relative contribution and quantitative contribution are different. Therefore, in this section we evaluated the quantitative contribution of each factor’s absolute effect, including AOD, on spring blooms using a PD plot. The contribution of predictors to fluctuating CHL was expressed as the CHL_<predictor>. The predictors included AOD (CHL_AOD), SST (CHL_SST), MLD (CHL_MLD), TP (CHL_TP), and CHL_c (CHLCHL_c). To quantitatively reveal the fluctuation of CHL according to the increase or decrease in AOD, we fixed the values of other predictors, except for AOD, during the entire period (i.e., an average of variables for March 2003 to 2021), as shown in Eq. (2). The fluctuating CHL according to the change in AOD was then determined using Eq. (1). For other predictors (e.g., SST, MLD, WD, TP, and CHL_c), the fluctuating CHL was simulated using a method similar to that of AOD.

Figure 6 shows the estimated CHL_RF through the PD plots. This analysis can explain the fluctuating CHL_RF responses to each variable. When AOD increased from 0 to 1, CHL_AOD increased from 1.31 mg m^-3 to 1.76 mg m^-3 (Figure 6A). When AOD was less than approximately 0.3, CHL_AOD slowly increased. However, the peak value of CHL_AOD appeared at 0.3 and the CHL fluctuation for AOD was approximately 1.40 mg m^-3. Hence, a low AOD during the COVID-19 lockdown period induced a slight variation in CHL_AOD. Therefore, through the PD plot for AOD, AOD reduction was presumed to cause weaker spring blooms in March 2020 than in other years. As a result, the AOD during the COVID-19 lockdown period (March 2020) recorded an abnormally low CHL (Section 3.1 and Figure 2). CHL_SST increased from 1.35 mg m^-3 to 2.20 mg m^-3 when SST changed from 5°C to 11°C (Figure 6B). From approximately 7°C to 8.5°C , where SST was low, phytoplankton grew slowly and even decreased in abundance (Xuan et al., 2011). However, SST gradually increased and reached 8°C to 12°C (the average SST in the CYS ranges from 8°C to 10°C ), showing a spring bloom peak (Lu et al., 2021). Moreover, the changing CHL pattern influenced by SST matched the generated spring bloom pattern conditions in the CYS (Xu et al., 2013; Zhao et al., 2019).

When MLD changed from 15 m to 85 m, CHL_MLD changed slightly from 1.57 mg m^-3 to 1.74 mg m^-3 (Figure 6C), thereby reflecting the mean MLD (~ 60 m) in the CYS. Furthermore, the CHL_MLD increased with MLD depth. In the CYS, nutrients that affect phytoplankton growth accumulate in the bottom water during autumn and winter. The accumulated nutrients are then moved from the lower layer to the upper layer by vertical mixing in the spring (Kim et al., 2000). Additionally, the YS is affected by anthropogenic aerosols in the atmosphere (Lin et al., 2005; Wei et al., 2015). However, comparing Figures 6A, C, the rate of increase in CHL by MLD was relatively weaker than that in the case of AOD. This result is supported by the data in Figure 2. Both 2020 and 2021 showed a trend of shallower MLD than the average; however, CHL showed the opposite trend (Section 3.2 and Figure 2) and AOD showed the same variability as CHL during the same period. Consequently, we could indirectly estimate that CHL was changed by AOD (Figures 2A, F). These results are consistent with those of previous studies. According to Moon et al. (2021), since 2010, the YS has been significantly affected by anthropogenic nitrogen deposition in the atmosphere. Further, Yoon et al. (2021) noted that most of the new nitrogen input in the YS was caused by atmospheric nitrogen deposition. Therefore, we found that the fluctuation in CHL was more affected by AOD than by MLD.

When WD ranged from − 1.09 × 10^-5 sec^-1 to 4.53 × 10^-6 sec^-1, CHL_WD fluctuated from 1.53 mg m^-3 to 1.68 mg m^-3 (Figure 6D). Negative divergence refers to the WD flow from the sea surface to the atmosphere. If the WD flow from the sea surface to the atmosphere decreases, WD flow from the atmosphere to the ocean (positive divergence) increases. As WD increased, CHL_WD decreased and then gradually increased (Figure 6D). The WD for March 2020 showed negative divergence. In addition, it revealed partially restricted industrial activities owing to the COVID-19 lockdown policy and showed that anthropogenic aerosol emissions decreased. As a result, the WD flow from the atmosphere to the sea surface weakened and the amount of aerosols deposited on the sea surface was quickly reduced. Conversely, in 2021, the flow from the atmosphere to the sea surface increased, which could have added abundant amounts of aerosols to the ocean. This process may act as a factor limiting the occurrence of spring blooms, thereby supporting the hypothesis that reducing AOD caused a decline in phytoplankton biomass. When TP increased from 0 mm to 0.15 mm, CHL_TP decreased from 2.42 mg m^-3 to 1.28 mg m^-3 (Figure 6E). These results suggest that phytoplankton growth decreased with precipitation. Furthermore, according to previous studies, phytoplankton growth is indirectly estimated to be limited due to the toxic effects of metal components, such as Cu and Pb, in precipitation (Liu et al., 2013). When CHL_c increased from 1.29 mg m^-3 to 3.55 mg m^-3, CHLCHL_C increased from 1.24 mg m^-3 to 4.07 mg m^-3 (Figure 6F). Because CHL_c represents the climatology of CHL in the CYS, CHLCHL_C increased with increasing CHL_c.

Based on the results of the quantitative contributions in Section 3.4, we created an equation that can express the CHL fluctuation for each predictor. Thus, this equation can estimate the amount of CHL fluctuation exhibiting different responses to each factor by summing the amount of change by the deviation of each factor in a year. Section 3.5 reveals the quantitative contribution of each predictor to CHL under natural conditions. The resulting equations are as follows:

where α_n is the empirically determined regression coefficient for each predictor (Figure 6; Table 3). Additionally, Δ CHL_<predictor> (difference CHL_<predictor>) represents the difference between x and x−1 , where x represents each year of the study periods 2021, 2020, 2019, … and 2004. Hence, Δ CHL_<predictor> was calculated for each variable variation simulated in the PD plot. Subsequently, the fluctuation CHL was calculated by summing the values of each Δ CHL_<predictor>.

Although the PD plot in Figure 6 explains the responsive function of each predictor to CHL changes, we cannot estimate specific contributions of each predictor in different years. Subsequently, the αn (empirical regression coefficients, 2aX + b, X is predictor) was defined by calculating the first derivation based on the empirical polynomial regression. The first derivative can represent the effect of a predictor on CHL (e.g., positive or negative linear trends) (Boelkins et al., 2018). All regression coefficients are listed in Table 3. We calculated the quantitative contribution of each predictor to CHL under natural conditions only in March 2020 (during the COVID-19 lockdown period). The α_n for each selected predictor was applied to Eq. (5). According to the variation in each predictor, the CHL variability was then calculated (Table 3). The difference in CHL with the variational predictor (Δ CHL_<predictor>) was calculated as the difference between the COVID-19 lockdown period and the average from March 2003 to 2021 (climatology). Δ CHL_AOD is the difference in CHL_AOD, which indicates the change between the COVID-19 lockdown period and climatology after fixing the variable conditions, except AOD, as in Section 3.4. The other variables were calculated using the same process. Using the α_n based on Eq. (5), Δ CHL_<predictor> was quantitatively calculated for each predictor. Figure 7 shows the effects of each variable on Δ CHL_<predictor> when it decreased by at least 90% and increased by 90% at most, compared to the climatology. The variation range for the predictor is the resulting value for all conditions in which each factor may change, including the COVID-19 lockdown period. Furthermore, Figure 7 shows the individual effect of the absolute contribution of each parameter to the CHL. As a result, Δ CHL_AOD decreased by 0.20 mg m^-3 due to an AOD less than that of the climatology. However, Δ CHL_SST increased by 0.22 mg m^-3 because SST was warmer than that measured in the same period in other years (Figure 7). Moreover, MLD became more shallow than that of the climatology, reducing Δ CHL_MLD by 0.01 mg m^-3. Additionally, the negative WD yielded a stable Δ CHL_WD, but less TP increased Δ CHL_TP by 0.02 mg m^-3. Similary, we estimated Δ CHLAOD in 2019 and 2021 using the linear model of the RF in Figure 7, showing as much as -0.05 mg m^-3 and + 0.13 mg m^-3, respectively. However, Δ CHLSST increased by 0.07 mg m^-3 and 0.18 mg m^-3 in 2019 and 2021, respectively. Thus, the linear model allows us to estimate the CHL chages due to the various predictors

From analyzing the individual effects of each factor, it was determined that the decreased AOD and shallower MLD, compared to those of climatology, affected the CHL reduction. Among them, the leading cause of the decreasing spring blooms during the COVID-19 lockdown period was the sharply reduced AOD from the shutdown policy in China. In addition, it was found that the increased SST and decreased TP effect on the increasing CHL. As a consequence, AOD and MLD were the causes that led to the decrease in CHL and, among them, the effect of AOD reduction was approximately 10 times greater than that of the shallower MLD. In addition, the quantitative (absolute) contribution of the individual factors to CHL, shown in Figure 7, was largest in SST, followed by those in AOD, MLD, TP, and WD. For comparison, they are listed in absolute value order without mathematical symbols. This differs from the relative contributions presented in Section 3.2, suggesting that the strength of spring blooms in the YS can be controlled mainly by SST and AOD. In addition to the relationship already revealed in Figure 7, the final results of this study, which is the quantitative relationship between CHL and each factor, which was limited in previous studies, are presented together.

Accordingly, we determined a causal relationship between AOD and CHL through time series analysis for each factor to investigate the cause of spring blooms that were weakened during the COVID-19 lockdown period. In addition, intense spring blooms were induced when additional AOD was supplied in March 2020 using the RF model. Finally, as a result of estimating the quantitative contribution of each factor to CHL, it was observed that the decrease in aerosols led to a decline in CHL. Overall, the weakened spring blooms during the COVID-19 lockdown period could have resulted from complex factors instigated by the temporary closure of industrial activities.