Independent components of spatial-temporal structure of chlorophyll a patterns in the upper layer of the north-western shelf of the Black Sea

In the present study, the results of independent component decomposition of satellitederived chlorophyll a (Chla) patterns for the north-western part of the Black Sea are presented. The study has been carried out on the basis of the DINEOF-reconstructed dataset of 8-day average log-transformed Chla (alChla) patterns for 1997-2016. The alChla patterns were decomposed into six independent components of its spatio-temporal variability in the north-western shelf of the Black Sea. The independent components reflect the spatial distribution of alChla anomalies which are likely to be formed under the influence of sea circulation factors driven by wind. The paper presents the results of the analysis of the intra-annual variability of independent components. The interpretation of the patterns of intra-annual independent components variability is given, taking into account the seasonal variability of the wind factor, the flow of the Danube, the Dnieper and Southern Bug rivers and the fact of modulation of independent components dynamics by seasonal phytoplankton succession. ‡ § © Kyrylenko N, Evstigneev V. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.


Introduction
The Black Sea is a semi-enclosed basin with limited water exchange through the Bosphorus Strait with the Mediterranean Sea. Strong salinity and density stratification of the Black Sea waters are due to the influx of fresh waters of the Danube, the Dnieper, the Dniester and the Southern Bug rivers on the one hand and salty Mediterranean waters coming into deep layers of the sea on the other hand. The main pycnocline, formed as a result of this, limits the ventilation of the lower layers of the Black Sea (below 150 m). Thus, only the upper oxygen-containing thin layer (60-200 m) provides the conditions for the existence of biota (Oguz et al. 2012, Pokazeev et al. 2020. These factors, along with a rather vast shelf zone in the north-western part of the Black Sea with its own characteristics of biogeochemical processes, determine the spatio-temporal heterogeneity of the functioning of the sea ecosystem (BEGUN et al. 2010) . Due to the shallow water and the influence of the large rivers ( Fig. 1), the circulation of water masses on the northwestern shelf has distinctive features compared to the deepwater part of the sea (Ivanov and Belokopytov 2011). Being an interesting study object, both from oceanological and biological points of view, this region of the Black Sea attracts the attention of many scientists. The recent studies have focused on drivers and seasonal variability of harmful hypoxia events (Capet et al. 2013).
One of the indicators of the marine ecosystem functioning is considered to be the content of the photosynthetic pigment -chlorophyll a (Chla) (McQuatters-Gollop et al. 2008, Nezlin North-western shelf of the Black Sea and the main rivers of the region. 2006). The concentration of Chla enables one to obtain data on the biomass of phytoplankton in sea surface layer by using the known ratio between Chla and the content of organic nitrogen in marine phytoplankton (Ragueneau et al. 2002). Data on Chla are also used to calculate the primary production of the oceans and seas and to validate the models of functioning of marine ecosystems in general, including in the Black Sea (Curran et al. 2018). Identification of the content of photosynthetic pigment by its fluorescence is a widely used method. Although this method is representative, it is quite labour-intensive, requiring carrying out field expeditions with the involvement of relevant specialists. Currently, existing satellite spectroradiometric scanners allow remote monitoring of the content and spatial distribution of Chla ( Suslin et al. 2018), based on the spectrum of upward radiation from the sea surface. The use of satellite information has become an integral part of the study of marine ecosystems and the factors conditioning their functioning (McQuatters-Gollop et al. 2008, Staneva et al. 2010, Curran et al. 2018. The aim of this study is to identify the main types of spatio-temporal structure of Chla patterns, which are formed under the influence of hydrodynamic factors of the marine environment in the surface layer of the north-western shelf zone of the Black Sea. The main source of information about Chla is satellite-sensing data, presented in different time resolutions. On the one hand, wind-induced circulation adjustment on the shelf is related to the synoptic timescale processes; therefore, satellite data of sub-monthly resolution should be used. On the other hand, the dataset of this temporal resolution may contain lost pixels (or grid cells) due to cloud cover, atmospheric correction algorithm errors etc. The array of 8-day Chla patterns seemed to us to be a reasonable compromise between these two factors determining the quality of the study results.
The paper is organised as follows. The data used in the study, as well as peculiarities of data transformation, are described in Section "Data". Section "Chla patterns decomposition technique" briefly describes a method for independent component analysis to be applied to Chla patterns. Section "Results and Discussion" presents the results of identifying independent spatio-temporal modes of Chla distribution for the north-western area of the Black Sea and joint analysis of the frequency of types and characteristics of wind circulation factors in a specified area.

Data
In order to carry out the study, we combine satellite data and grid data of re-analysis of wind patterns over the north-western part of the Black Sea. The time period to be analysed was limited by the span of available satellite data for this region.

Satellite Data
In this work, a dataset of the 8-day average Chla fields for the north-western region of the Black Sea was used. The dataset is presented in a single spatial resolution grid of 4.63 km and encompasses the period 1997-2016. These data have been produced in the framework of the GlobColour project of the European Space Agency (ESA) using the datasets of SeaWiFS, MODIS, MERIS and VIIRS scanners. Chla patterns from different scanners were combined by the Garver-Siegel-Maritorena (GSM) method (Maritorena and Siegel 2005) and can be downloaded from www.globcolour.info. Satellite-derived Chla is obtained from the ocean colour in different spectral channels using special algorithms, depending on the type of surface waters (Suslin et al. 2018). The water of the Black Sea corresponds to the type (case-2) where the surface colour depends on Chla, as well as on concentrations of dissolved and suspended substances. The use of the data generated by case-1 water recalculation algorithm leads to an overestimation of the north-western Black Sea Chla patterns in absolute value because this algorithm is applicable to water, the colour of which is determined only by the Chla content. However, the aim of this work is not to compare the absolute values of the Chla content with in situ measurements. As in Nezlin (2006), we used satellite sounding data as an indicator of the spatio-temporal variability of Chla in the studied region. Considering the availability of case-1 data for the longest period of time (since 1997), their use allows us to more reliably establish the modes of Chla distribution.
The 8-day Chla arrays, as a rule, contain lost pixels due to a number of factors (cloud cover, errors of the atmospheric correction algorithm etc.). The gaps were statistically reconstructed using empirical orthogonal functions according to the DINEOF method (Alvera-Azcárate et al. 2005), which makes it possible to restore integrity of the dataset with a significant number of data gaps.
The obtained Chla dataset was preprocessed to increase the reliability of the analysis results. The pigment content in the sea surface layer is a substantially positive log-normally distributed random variable, not only in the Black Sea (Nezlin 2006), but also in other seas and oceans (Wang and Liu 2013). In this regard, logarithmically transformed (natural logarithm) fields of Chla were used. Moreover, Chla variability may differ by several orders of magnitude for the estuaries of the region's main rivers -the Danube River, the Dnieper, the Dniester River -and in the areas adjacent to the deep-sea areas of the Black Sea (Staneva et al. 2010, Kyrylenko and Evstigneev 2017. To minimise the influence of the heterogeneity of Chla distribution, the initial data were converted into absolute anomaly fields where is the average long-term value of the log-transformed Chla value in each pixel of the field. The alChla data are grouped in the form of an observation matrix X , where m = 1, ..., M is the spatial index (pixel index) and t = 1, ..., T is the time index.

Re-analysis Data
In this paper, to characterise the external environmental factors and their variability, a daily data array of re-analyses ERA-Interim (at grid nodes of 0.25º x 0.25º) consisting of wind components (U, V), tangential wind stress in zonal (τ ) and meridional directions (τ ) was t×m X Y involved. In addition, to characterise the wind circulation of water in the surface layer of the Black Sea, the wind stress curl was calculated from the data of τ , τ .
Eight-day fields, corresponding to Chla data, were calculated from the re-analysis data by simple averaging.

Chla patterns decomposition technique
Analysis of multivariate time series, as a rule, is carried out to reveal hidden regularities associated with factors impacting the system under study and manifesting themselves as modes both in time and in space (Hannachi et al. 2007). There are many techniques for dimension reduction and decomposition of multivariate data arrays, but the principal component analysis (PCA), singular value decomposition (SVD) and rotational techniques are considered to be the most popular in biological and geophysical sciences. All these methods are based on several assumptions regarding the properties of a multivariate random variable to be decomposed: -the resulting space-time array is a linear mixture of principal components corresponding to factors impacting the system under study; -the distribution of stochastic variables, associated with principal components, are normally distributed random variables.
In fact, alChla fields have non-Gaussian distribution. Fig. 2 shows the maps of the spatial distribution of high-order statistics -a skewness and a kurtosis. Despite the logarithmic transformation of the Chla values, the kurtosis is significantly larger than 0 for the entire domain, i.e. alChla distribution is peaked. For the skewness, the results are not so unambiguous. The insignificant asymmetry of alChla distribution is characteristic for the central region of the north-western shelf of the Black Sea and to the south of the Danube Delta. In the vast coastal zones and the regions adjacent to the deep-water area of the sea, the skewness is significantly less than 0. Thus, the calculation proves in favour of a non-Gaussian distribution of the random variable of alChla. Moreover, factors impacting Chla patterns may also be non-Gaussian in nature.
Consequently, application of standard techniques for decomposition of multivariate Chla dataset is incorrect. In this work, the recently-developed method of decomposition into independent components is used for these purposes. Independent component analysis (ICA) (Aires et al. 2000) is based on assumptions about the non-Gaussian nature of the extracted independent components (ICs) and their statistical independence. ICA does not require the additivity and orthogonality of the components -the condition that can more accurately reflect the actual external and internal factors affecting the Chla distribution in shelf zone.

X Y
The idea of ICA decomposition is to treat observation matrix X as a mixture of N ≤ M independent "sources" of signals . The requirement for component independence is more stringent than their non-correlation (as one assumes when applying PCA). For independent random variables , the joint probability density function is expressed as a product of univariate functions. In case of non-Gaussian components, uncorrelatedness of s is not obvious. That is why any estimate of deviation of empirical probability density function from the Gaussian one may be considered as a measure of independence in ICA. According to the FastICA algorithm to be applied in this study, negentropy is considered to be a measure of independence (Hyvärinen 1999).
The independent components to be extracted represent a certain type of evolution of the alChla patterns or the so-called "prototype" of behaviour (Richman 1981). Interpretation of independent modes in terms of environmental factors should be based on an analysis of temporal variability of independent components and the geographic localisation of its spatial maxima/minima. In the simplest case, one factor (physical or biological) determines the specific dynamics of the system under study which manifests itself in different geographical areas of the domain. A more complex case implies that dynamic factors can have different modes manifesting themselves in different geographical areas. In this case, a number of behaviour "prototypes" may be associated with the same phenomenon. We consider several physical phenomena (wind over domain and Danube+Dnieper rivers inflow), which are believed to affect the dynamics of pigment content in the surface layer of the Black Sea north-western shelf.

Results and Discussion
Six independent components of alChla patterns were determined by the ICA method applied to a dataset of 8-day average satellite data on the content of Chla in the surface layer of the north-western shelf of the Black Sea. The choice of the number of types is ambiguous as there is no universal way to determine their optimal number. However, we focused on the results of Chla patterns decomposition into empirical orthogonal functions (EOF). Calculations shows that up to six empirical orthogonal modes account for more than 80% of the field total variance. The variance of the rest of the EOFs has values less than 2% and their loading maps are much distorted and badly distinguished. A comparison of the IC scores with their physical factors temporal variability allowed us to give an interpretation of the possible reasons for their formation.  (Sobotka et al. 2021). Fig. 3b shows the annual distribution of mean annual values of the IC1 with a clear positive peak in May-June and a minimum in late autumn and winter. An l×m i obvious factor determining such alChla dynamics is river run-off (especially the Danube River), albeit the peak of pigment content is shifted in relation to the run-off maximum (spring flood) by one month (Kyrylenko and Evstigneev 2017). Fig. 3c illustrates the joint seasonal variation of the average IC1 values with the normalised deviation of flow volume of the Danube River averaged over the period 1997-2016. Obviously, these curves coincide with a lag-1 month at least for the warm half of the year.

Independent component IC1
The peak of flood for the rivers of the north-western shelf, in particular the Danube and Dnieper Rivers, falls on April-May. With an increasing volume of fresh river waters carried to the shelf, the influx of nutrients assimilated by the phytoplankton within 1-2 months increases. The spring-summer peak of Chla is well-recognised by many studies (McQuatters-Gollop et al. 2008, Yunev 2011, Finenko et al. 2014, Kyrylenko and Evstigneev 2017, being accomplished using in situ measurements (Demidov 2008), as well as Chla derived from the satellite data by means of a regional algorithm for retrieval pigment concentration in the surface layer from ocean remote-sensing reflectance (Suslin and Churilova 2016). The shift between river inflow and the Chla peak can be explained by seasonal evolution of vertical stratification of the water column and phytoplankton succession. According to Capet et al. (2013), increased river discharge in spring leads to haline stratification of shelf waters. The subsequent heating of the surface waters enhances stratification and, by May, the mixed layer depth takes its minimum values. Wellwarmed and nutrient-rich waters favour development of large-sized species of dinoflagellates becoming dominant subsequent to diatoms prevailing in the previous winter-spring cold period (Mikaelyan et al. 2018).
Despite all the variety of factors affecting Chla during this period, we will focus on wind circulation over the area under study. Maximum correlation between 8-day mean values a b In the autumn-winter period, there is a mismatch between the course of IC1 and the run-off volume. In our opinion, this is due to the dominant influence of the wind factor. Fig. 3c depicts normalised wind curl anomalies over the north-western region according to ERA Interim data for the period of 1997-2016. It can be seen that, in the spring-summer period, the anticyclonic wind vorticity increases over the region which contributes to the development of a similar sea-surface circulation. As a result, this leads to the expansion of the spread area of the fresh-water to the area of the north-western shelf. In the period from October to December, the vorticity of the wind field becomes positive and determines the cyclonic circulation of seawater on the shelf. As a result, surface water from the deep-sea area of the Black Sea, enriched with Chla, flows into the shelf . IC1 describes about 26% of the total variance of the Chla dataset or its log-anomalies alChla to be more precise.

Independent component IC2
IC2 describes about 13.5% of the total variance of alChla log fields. The analysis of the base temporal function of IC2 shows an obvious intra-annual variability -the presence of two peaks in winter (February-March) and in autumn (September-November). On the loading map, these peaks correspond to two unrelated areas of localisation of positive values of the IC2 decomposition coefficients. The first narrow area stretched along the coastline near the mouths of the main rivers (Fig. 4). The second area with maximum positive a values is located on the slope of the continental shelf in the south-eastern part of the domain.
In order to determine hydrodynamical processes that form this type of a distribution, it is probably worth noting the wind factor under the influence of which sea surface currents develop both in the coastal region and on the slope of the north-western shelf of the Black Sea. This assumption is confirmed by correlation estimates (see alsoFig. 4c). The same as in Fig. 3, but for IC2.
formation of negative anomalies in the content of Chla in the shelf zone, while in the deepsea region, a winter maximum of pigment content can be observed (Finenko et al. 2014).
In the areas of river water spreading at depths of less than 50 m, winds of 4 m/s and more create a steady wind current in the upper layer and a countercurrent in the lower horizons. The boundary between these two layers is a thermocline layer with a sharp change in temperature and salinity. With an increase in wind speed up to 8 m/s, the wind current covers the entire column of water from the surface to the bottom. Long-lasting and directionally stable winds with a speed of 13-14 m/s create a single-layer flow with a speed of 1 m/s in the surface layer and down to 0.2-0.3 m/s at the bottom (at depth of 10-20 m).
During autumn-winter, convection surface water mixes with the bottom water, resulting in a gradual cooling. Ongoing convection until late autumn causes surface water to be completely transformed into bottom water over the entire sea area under study. In the estuaries, river water continues to flow into the sea providing an uninterrupted input of nutrients, mixing with and transforming into bottom water within the estuaries or in the immediate vicinity. It is worth noting that, according to the expedition data, a vertical distribution of nutrients has specific regional features (Oguz et al. 2002). For example, in the plume region, an increased content of phosphates and nitrates is observed only in a thin upper layer (Kukushkin and Parkhomenko 2018).
The structure of the IC2 pattern corresponding to the autumn peak (September -November) (Fig. 4b) is largely determined by phytoplankton succession. In autumn, diatoms (Sceletonema costatum, Chaetoceros wighamii), green algae and cyanobacteria are known to develop massively in estuarine areas (Finenko et al. 2008). The development of cyanobacteria is driven by low salinity which can favour their blooming. Positive a loadings in the shelf/deep-sea boundary region over the continental slope (isobath 200 m) (see Fig. 4a) may be associated with propagation of deep-sea water enriched with Chla due to the autumn secondary peak (Chu et al. 2005).
Intra-annual variability of IC2 is characterised by the prevalence of its negative values in summer. Thermohaline stratification of the sea water takes place in summer significantly complicating the upward transport of nutrients into the upper layer above the thermocline. Lack of nutrients leads to phytoplankton diversity and biomass decline in the upper layer (McQuatters-Gollop et al. 2008, Demidov 2008, Finenko et al. 2014.

Independent components IC3, 4 and 5
The IC 3-5 components describe about 23% (9.3%, 8.2% and 5.4%, respectively) of the total variance of alChla patterns. All three ICs are depicted in Figs 5, 6, 7 since the meridional wind component is the main factor driving such "behaviour prototypes" in our opinion (see Table 1). The shallow nature of the shelf zone contributes to the rapid reaction of the current structure to wind stress (Ivanov and Belokopytov 2011). The maximum mean IC scores (Figs 5b,6b,7b) in the spring months are primarily associated with the influence of the total river inflow. The features of distribution of river waters on the shelf under the influence of drift currents determined the spatial distribution of a on the loading maps (see Figs 5c, 6c, 7c).  Table 1.
Maximum absolute values of the correlation coefficient between wind characteristics and independent component ICs. Coefficients are significant at 5% level.

Figure 5.
The same as in Fig. 3, but for IC3.
The vorticity of the wind field in the spring months is insignificant on average, whereas the meridional wind component reaches its maximum. In particular, under the influence of the meridional component of the wind, the Danube waters can spread to the shelf north of Snake Island up to the Dniester Estuary (see Fig. 5a). On the other hand, in early spring, the wind (especially its meridional component) contributes to the formation of a Danube anticyclone, which extends to the east by summer, spreading along the entire shelf and providing anticyclonic circulation of waters enriched with nutrients (see Figs 6a, 7a). That might happen due to baroclinic instability caused by intensive desalination, but many experts tend to unequivocally believe in the wind origin of this eddy (Ivanov and Belokopytov 2011).
In summer, the wind curl takes negative values which favours anticyclonic circulation in the north-western shelf zone. Thus, loading maps in the summer period reflect different variants of evolution of alChla patterns under the influence of the anticyclonic vorticity of the wind field ).
It should be noted that mesoscale variability in the north-western shelf plays an important role in the structure of "prototypes" IC3, 4 and 5. As a result of interaction of the Black Sea a b c Figure 6.
The same as in Fig. 3, but for IC4.
Rim Current and Sevastopol Anticyclone which form a system of currents on the outer boundary of the shelf, the open sea waters expand to the north-western shelf zone. In addition, it was found out that, with considerable variability in surface currents under the influence of the rivers' freshwater inflow from the Dnieper-Bug Liman, an anticyclonic eddy was also formed near Odessa. Such a process development pattern is quite likely. For example, in ), a possible development of the Danube river plume, its transferring and confluence with the Rim Current, was studied. The authors showed that the recurrence of anticyclonic eddy near Odessa is about 15%.

Independent component IC6
The IC6 component describes 6.7% of variance of alChla patterns. Despite the insignificant contribution, IC6 cannot be ignored since it reflects the dynamics of changes in the Chla pattern on the north-western shelf in the winter-spring period.
Large positive values of the wind curl indicate an increase in cyclonic circulation in the region. This type of circulation in the cold half-year determines the intensification of the local anticyclonic eddies. In the western part of the shelf, the Danube anticyclonic eddies are formed in early spring, which contribute to the spread of river waters. The region of its formation coincides with the western area of maximum matching of loadings a on the map (see Fig. 8a). In the eastern part, during cyclonic circulation on the shelf, the Kalamitsky a b c j Figure 7.
The same as in Fig. 3, but for IC5. and Karkinitsky anticyclones intensify, contributing to non-propagation of water into the central part of the shelf. In addition, the increased pigment content in these regions is associated with the influx of pigment-enriched waters from the open part of the Black Sea due to the interaction of the Rim Current and the Sevastopol anticyclone (Finenko et al. 2014).
In summer, the IC6 loading map has an opposite distribution of maxima and minima (negative IC6 values in Fig. 8b) which is more likely to reflect the influence of continental waters from the Dnieper-Bug Estuary and anthropogenic impact on phytoplankton development. It is known that increased content of Chla and phytoplankton biomass occur when the plume of transformed river freshwater enters the Odessa region from the Dnieper-Bug Estuary in spring. High concentration of nutrients in the river water contributes to the primary production of the organics by phytoplankton. This organic substance is deposited in sediment layers resulted in regeneration of nutrients and eutrophication of the water in this area in the summer season.

Conclusions
This study has been carried out on the basis of reconstructed arrays of Chla fields of 8-day averaging (using the DINEOF method) during the period of 1997-2016. The Chla fields were decomposed into six independent components of spatio-temporal organisation of a b c Figure 8.
The same as in Fig. 3, but for IC6.