1. INTRODUCTION / Uvod
Primary production (PP) in the Adriatic Sea is spatially heterogeneous. The northern Adriatic is shallow, strongly influenced by riverine inputs and generally characterized by higher productivity, whereas the deeper central and southern Adriatic are more oligotrophic and characterized by lower production [1-6].
A major methodological step in the study of marine PP was the introduction of the radiocarbon ¹⁴C method by Steemann Nielsen in 1952. This method became one of the most widely used approaches for measuring marine primary production and remains the most frequently applied method in Adriatic PP studies. Soon after its introduction, the method was implemented by scientists from the Institute of Oceanography and Fisheries in Split. Regular in situ PP measurements in the middle Adriatic started in April 1962 at the open-sea station Stončica, near the island of Vis, and at the coastal station Kaštela Bay, establishing the longest continuous ¹⁴C-based primary production time series in marine science [7-10].
The importance of Stončica is linked not only to the length of the record, but also to its position in the middle Adriatic. This station represents open-sea conditions, but it is still part of a basin affected by water-mass exchange, circulation changes and occasional external inputs [6]. Long-term observations have shown that PP at Stončica did not follow a simple monotonic trend. Instead, the record includes periods with different productivity regimes [11], including phases of increased and reduced PP [8,12-15]. These changes have been related to anthropogenic pressure [8,16], atmospheric variability [17,18], Mediterranean water ingression [6,19] and the Adriatic-Ionian Bimodal Oscillating System (BiOS) [15,20-22]. Previous long-term studies in the middle Adriatic have further shown that variability in PP should be considered within a broader hydrographic and climatic context. Changes in water-mass properties, atmospheric forcing and biological and chemical conditions may operate across different temporal scales, linking local productivity variability with wider Adriatic and Adriatic–Ionian processes [11,14,17–20]. This broader context is particularly relevant when interpreting long-term PP variability across areas that differ in their exposure to coastal, riverine and open-sea influences.
The scientific value of the Stončica time series lies primarily in its exceptional duration. Measurements were conducted repeatedly throughout the year, often at an approximately monthly scale, although the exact number and timing of observations varied among years. In some years the sampling was close to monthly, whereas in others it was closer to seasonal coverage [10]. This means that the main limitation of the record is not its length, but its temporal irregularity. Events shorter than the sampling interval, including short phytoplankton blooms, wind-driven mixing episodes, episodic atmospheric deposition or mesoscale processes, could be missed or only partly represented. In addition, in situ PP measurements are spatially limited and cannot resolve gradients between coastal, river-influenced and more open-sea areas across the central and southern eastern Adriatic.
Satellite ocean-colour observations provide a useful complementary approach to these limitations. Ocean-colour sensors can observe large marine areas repeatedly, although the effective availability of usable observations is reduced mainly by cloud cover and quality-control constraints. For long-term analyses, merged multi-sensor products such as ESA OC-CCI are particularly useful because they combine observations from several satellite missions into a harmonised record [23]. In this study, gap-free monthly satellite-derived PP estimates from OC-CCI Version 4.2 were used for the period 1998–2020. These values should not be interpreted as single measurements equivalent to one in situ ¹⁴C incubation. They are monthly composite estimates derived from satellite observations and bio-optical modelling. Satellite-derived PP therefore provides regular monthly temporal coverage and broad spatial coverage, but it remains an indirect, algorithm-dependent estimate of carbon fixation.
Although long-term in situ observations have provided valuable information on temporal changes in PP in the Adriatic, spatially consistent comparisons of long-term PP variability among contrasting coastal, river-influenced and offshore areas of the central and southern eastern Adriatic remain limited. By analysing five areas within the same satellite-derived dataset, this study provides a common framework for distinguishing shared regional variability from differences among areas exposed to different environmental settings.
The aim of this study was to assess spatial, seasonal and interannual variability of satellite-derived PP in five selected areas of the central and southern eastern Adriatic Sea: Split Channel, Stončica, Neretva, Dubrovnik and Albania. These areas include urban coastal, river-influenced, southern coastal and more open-sea conditions. Monthly and annual PP values were analysed using spatial comparisons, seasonal cycles, linear trend analysis and Dynamic Factor Analysis (DFA). The study examines whether satellite-derived PP can complement the historical ¹⁴C record by providing a regular monthly framework for detecting common regional variability and area-specific responses across the analysed areas.
2. MATERIALS AND METHODS / Materijali i metode
2.1. Study area and sampling sites / Područje istraživanja i točke uzorkovanja
This study was conducted in selected areas of the central and southern eastern region of the Adriatic Sea during the period from 1998 to 2020. Five study areas were chosen for analysis, namely Split Channel, Stončica, Neretva, Dubrovnik, and Albania (Figure 1). The study was restricted to selected areas of the central and southern eastern Adriatic. The northern Adriatic was not included because it represents a distinct shallow and strongly river-influenced subsystem, characterized by substantially different hydrographic and productivity conditions from the deeper central and southern Adriatic [5,6]. The present study was therefore designed to compare coastal, river-influenced and offshore conditions within the central and southern part of the eastern Adriatic rather than to represent the entire basin. These study areas have been chosen to enable comparison of satellite-derived productivity patterns within the areas impacted by diverse environmental conditions, including riverine, coastal and open-sea conditions. Split Channel is an enclosed urban coastal area, where the PP could be influenced by coastal, anthropogenic and water exchange factors. Stončica is an offshore study area in the middle Adriatic Sea that is less affected by the terrestrial environment, oligotrophic in nature. Neretva is a shallow coastal site that is influenced by the river Neretva, which would influence the PP through nutrients and freshwater inputs. Dubrovnik is a coastal study area in the southern part of the Adriatic Sea, while Albania is a coastal site in the southeastern part of the Adriatic Sea, which is influenced by the river Bojana/Buna. PP in the central and southern Adriatic generally follows a pronounced seasonal cycle controlled by changes in light availability, vertical mixing, nutrient supply and water-column stratification. Winter mixing replenishes nutrients in the upper water column, while increasing irradiance and progressive stabilization of the water column favour enhanced phytoplankton production during spring. During summer, strong thermal stratification limits the upward supply of nutrients to the euphotic layer, resulting in lower PP. The magnitude and timing of this seasonal cycle may differ between offshore and coastal areas, particularly where freshwater and nutrient inputs modify local conditions [1,4,5].

Figure 1 Location of the five study areas in the central and southern eastern Adriatic Sea: Split Channel, Stončica, Neretva, Dubrovnik, and Albania. Coloured rectangles indicate the spatial extraction boxes used to obtain satellite-derived primary production data. The inset shows the position and extent of the study area within the Adriatic Sea.
Slika 1. Lokacija pet istraživanih područja u srednjem i južnom istočnom Jadranskom moru: Splitski kanal, Stončica, Neretva, Dubrovnik i Albanija. Obojeni pravokutnici predstavljaju prostorne okvire unutar kojih su prikupljeni satelitski podaci o primarnoj proizvodnji. Umetnuti prikaz pokazuje smještaj i opseg istraživanog područja u Jadranskom moru.
Source: Authors’ own work using Natural Earth data
2.2. Data source and data collection / Izvor i prikupljanje podataka
Satellite-derived PP data were obtained from the BICEP/NCEO global marine phytoplankton primary production dataset, derived from ESA Ocean Colour Climate Change Initiative (OC-CCI) Version 4.2 monthly products, at a horizontal spatial resolution of 9 km, downloaded from https://www.oceancolour.org [24]. PP was estimated using a spectrally resolved primary-production model that incorporates the vertical structure of phytoplankton and represents photosynthesis as a function of irradiance using a two-parameter photosynthesis–irradiance relationship. We did not calculate PP from the underlying ocean-colour variables or apply any additional regional calibration for the Adriatic Sea but used the precomputed satellite-derived PP product.
The analysed dataset consisted of monthly mean PP values for the period from January 1998 to December 2020, expressed as depth-integrated daily production in mg C m⁻² d⁻¹. The underlying Level-3 ocean-colour data are not gap-filled. Therefore, PP was analysed over predefined spatial extraction boxes rather than at individual pixels. For each study area, all available valid PP grid cells within the minimum and maximum longitude and latitude limits specified in Table 1 were used to calculate the monthly spatial mean, without imposing a minimum number of valid pixels. Grid cells flagged as invalid or containing missing PP values were excluded from calculation, and the monthly spatial mean was calculated only from the valid PP grid cells available within each extraction box. No temporal interpolation or additional gap filling was applied. Consequently, the final dataset contained one monthly mean PP value for each study area from January 1998 to December 2020.
The use of spatial extraction boxes rather than individual pixels reduces sensitivity to the spatial resolution of the satellite product and avoids overinterpreting variability associated with individual pixel locations. Small-scale local differences, such as those within the Malostonski Bay area, cannot be reliably resolved from the available satellite data; therefore, these areas were interpreted as parts of broader regional systems, in this case the Neretva-influenced coastal area. Satellite-derived PP estimates in coastal waters should nevertheless be interpreted with caution. The 9-km spatial resolution may smooth small-scale coastal gradients, while optically complex conditions associated with suspended material, coloured dissolved organic matter and terrestrial inputs may introduce additional uncertainty in ocean-colour-based PP estimates. The extracted values were subsequently organized by study area, year and month.
Table 1 Geographic boundaries of the spatial extraction boxes used for satellite-derived primary production data. Longitude and latitude ranges represent the minimum and maximum coordinates of each extraction box for the five analysed areas in the central and southern eastern Adriatic Sea: Split Channel, Stončica, Neretva, Dubrovnik, and Albania.
Tablica 1. Geografske granice prostornih okvira korištenih za satelitske podatke o primarnoj proizvodnji. Rasponi geografske dužine i širine predstavljaju minimalne i maksimalne koordinate svakog okvira za izdvajanje za pet analiziranih područja u srednjem i južnom istočnom Jadranskom moru: Splitski kanal, Stončica, Neretva, Dubrovnik i Albanija.
The analysis was restricted to 1998–2020 because this corresponds to the full temporal coverage of the 9-km primary-production product used in this study. More recent OC-CCI ocean-colour observations are available, but an equivalent continuation of the same primary-production product beyond 2020 was not available.
2.3. Data processing and statistical analysis / Obrada podataka i statistička analiza
Monthly satellite-derived PP data were analysed for each study area to assess spatial, seasonal and interannual variability. Data processing, statistical analyses and visualization were performed in R. After extracting the corresponding year and month from the satellite time series, monthly and annual mean PP values were calculated using the dplyr and lubridate R packages. Graphical outputs were produced with ggplot2, spatial visualization was prepared using sf and rnaturalearth, and DFA [25] was performed using the MARSS package [26].
Spatial variability was assessed from the distribution of monthly PP values at each study area and presented using boxplots. Seasonal variability was analysed by grouping monthly PP values by calendar month and calculating mean monthly PP for the entire 1998–2020 period. For statistical comparisons of PP among study areas and seasons, Kruskal–Wallis tests were applied to the monthly PP values. Seasons were defined as winter (December–February), spring (March–May), summer (June–August) and autumn (September–November). When a significant Kruskal–Wallis result was obtained, pairwise Wilcoxon rank-sum tests were used for post-hoc comparisons, with p-values adjusted using the Benjamini–Hochberg procedure. Differences in seasonal patterns among study areas were further examined using a two-way linear model including study area, season and their interaction (PP ~ study area × season). The significance of the study area × season interaction was evaluated using analysis of variance.
Interannual variability was evaluated using annual mean PP values calculated separately for each study area. Standardized annual anomalies (z-scores) were calculated separately for each study area as
,
where x is the annual mean PP, x̄ is the area-specific mean PP for the 1998-2020 period, and SD is the corresponding standard deviation. The resulting standardized anomalies are dimensionless. Linear trends were estimated independently for each study area using linear regression in R, with year as the explanatory variable and annual mean PP as the response variable. Trend coefficients, R² values and p-values were used to describe the direction, explanatory power and statistical significance of the trends, with p < 0.05 considered statistically significant. Because the trend analyses were based on annual time series, observations from successive years may not be fully independent. Temporal autocorrelation was not explicitly modelled in the simple linear regressions; consequently, the associated p-values were interpreted cautiously and the trends were considered primarily as descriptions of long-term directional change. Separate trend models were fitted for each of the five study areas, and the resulting p-values therefore represent area-specific tests rather than a single overall test of temporal change. No adjustment for multiple testing was applied across the five separate area-specific trend models; therefore, individual trend p-values should be interpreted with appropriate caution.
Area-specific long-term changes were further visualized by plotting annual mean PP values together with fitted linear regression lines. Finally, DFA was applied to the annual PP time series to identify temporal variability shared among the analysed study areas. DFA is a multivariate state-space method in which multiple observed time series are represented by a smaller number of unobserved common trends and area-specific variability [25]. In the present study, DFA was used to extract the common temporal pattern shared among the five study areas and to distinguish this regional signal from area-specific behaviour. The analysis was performed using the MARSS package in R [26].
3. RESULTS / Rezultati
Spatial and seasonal variability of satellite-derived PP is presented in Figure 2. Across all monthly values, the highest mean PP was recorded at Albania (450 mg C m⁻² d⁻¹), followed by Split Channel (422 mg C m⁻² d⁻¹) and Neretva (414 mg C m⁻² d⁻¹). Lower mean values were obtained for Dubrovnik (350 mg C m⁻² d⁻¹) and Stončica (328 mg C m⁻² d⁻¹). The relative variability was highest at Neretva, with a coefficient of variation (CV) of 31.5%. Albania and Dubrovnik had similar CV values, 28.6% and 28.4%, respectively, whereas Split Channel and Stončica showed lower relative variability, with CV values of 26.5% and 26.8%.
The seasonal cycle was similar at all study areas, with a spring maximum and a summer minimum (Figure 2b). Maximum PP occurred in May at all study areas except Neretva, where the maximum was recorded in April. The highest seasonal maximum was observed at Albania in May (679 mg C m⁻² d⁻¹), followed by Split Channel in May (624 mg C m⁻² d⁻¹), Dubrovnik in May (538 mg C m⁻² d⁻¹), Neretva in April (532 mg C m⁻² d⁻¹) and Stončica in May (487 mg C m⁻² d⁻¹). Minimum PP occurred in August at all study areas and ranged from 223 mg C m⁻² d⁻¹ at Stončica to 272 mg C m⁻² d⁻¹ at Albania. Following the August minimum, PP increased again in September at all analysed study areas, producing a secondary late-summer/early-autumn increase in the seasonal cycle. The seasonal amplitude was highest at Albania (407 mg C m⁻² d⁻¹) and lowest at Stončica (264 mg C m⁻² d⁻¹).
The Kruskal-Wallis test showed significant differences in PP among study areas (χ² = 204.94, df = 4, p < 0.001) and among seasons (χ² = 698.48, df = 3, p < 0.001). Pairwise post-hoc comparisons showed that most study area pairs differed significantly, except Neretva and Split Channel. Seasonal post-hoc comparisons showed significant differences among most seasonal pairs, except autumn and winter. The study area × season interaction was also significant in the two-way linear model (F = 7.82, p < 0.001), showing that the magnitude of the seasonal cycle differed among study areas.

Figure 2 Spatial and seasonal variability of satellite-derived primary production at the five study areas in the central and southern eastern Adriatic Sea: Split Channel, Stončica, Neretva, Dubrovnik and Albania, during the period 1998–2020. (a) Distribution of monthly satellite-derived PP values in the five analysed areas during 1998–2020; (b) mean seasonal cycle calculated from monthly PP values grouped by calendar month over the 1998–2020 period. Values are expressed in mg C m⁻² d⁻¹.
Slika 2. Prostorna i sezonska varijabilnost primarne proizvodnje temeljene na satelitskim podacima na pet istraživanih područja u srednjem i južnom istočnom Jadranu: Splitski kanal, Stončica, Neretva, Dubrovnik i Albanija, tijekom razdoblja 1998. – 2020. (a) Raspodjela mjesečnih vrijednosti PP-a temeljenih na satelitskim podacima u pet analiziranih područja tijekom 1998. – 2020.; (b) prosječni sezonski ciklus izračunat iz mjesečnih vrijednosti PP-a grupiranih prema kalendarskom mjesecu u razdoblju 1998. – 2020. Vrijednosti su izražene u mg C m⁻² d⁻¹.
Source: Authors’ own work based on BICEP/NCEO data
Interannual variability of annual mean PP is shown in Figure 3. Annual PP values fluctuated throughout the 1998–2020 period, but the magnitude and timing of these fluctuations differed among study areas. Albania generally had the highest annual PP values, while Stončica consistently showed the lowest values. Several years were characterized by pronounced standardized anomalies relative to the area-specific long-term mean. The strongest positive standardized anomaly was recorded at Albania in 2010 (z = +3.20), while positive anomalies were also observed in the same year at Split Channel (z = +2.07) and Neretva (z = +1.64). In contrast, reduced PP was particularly evident in 2018 at Neretva (z = −2.26) and Split Channel (z = −1.73).
The study areas did not show identical interannual behaviour. Stončica differed from the other study areas, with its strongest negative anomaly in 2001 (z = −1.66) and its strongest positive anomaly in 2020 (z = +2.31). Dubrovnik showed a more moderate range of interannual variability, with the highest positive anomaly in 1998 (z = +1.83) and the lowest negative anomaly in 2017 (z = −1.67). These differences indicate that, although some years were characterized by a common signal across several study areas, local or area-specific variability remained important.
Correlation analysis of annual mean PP supported this pattern. All correlations among study areas were positive, indicating that the study areas tended to vary in the same general direction over time. The strongest correlation was found between Neretva and Split Channel (r = 0.81), followed by Albania and Dubrovnik (r = 0.77). In contrast, Stončica showed weaker correlations with Neretva (r = 0.20) and Split Channel (r = 0.22), suggesting a more distinct interannual behaviour at the more open-sea study area. These results indicate partial synchronization of PP variability among study areas but also justify the subsequent use of DFA to separate the common regional signal from area-specific variability.

Figure 3 Interannual variability of satellite-derived primary production at the five study areas in the central and southern eastern Adriatic Sea during the period 1998–2020. Values are expressed as annual mean primary production in mg C m⁻² d⁻¹. Arrows indicate the year of the strongest standardized anomaly at each study area, with the corresponding z-score shown next to the label.
Slika 3. Međugodišnja varijabilnost primarne proizvodnje temeljene na satelitskim podacima na pet istraživanih područja u srednjem i južnom istočnom Jadranu tijekom razdoblja 1998. – 2020. Vrijednosti su izražene kao prosječna godišnja primarna proizvodnja u mg C m⁻² d⁻¹. Strelice označavaju godinu s najjačom standardiziranom anomalijom na svakom istraživanom području, a odgovarajući z-score prikazan je uz oznaku.
Source: Authors’ own work based on BICEP/NCEO data
DFA was applied to the annual PP time series to identify the common interannual signal shared among the analysed study areas. The extracted common trend showed pronounced interannual variability during the 1998–2020 period (Figure 4, Table 2). Positive DFA values indicate years in which the common PP signal was above its long-term average, whereas negative values indicate below-average common PP conditions.
The strongest positive common signal was observed in 2010, with a DFA value of 1.75. Positive phases were also recorded in 1998 (1.21) and 2004 (0.87), indicating years with above-average common PP conditions. In contrast, the lowest common signal occurred in 2017 (−1.03), followed by 2018 (−0.86), indicating a period of reduced common PP conditions. The common trend did not show a monotonic increase or decrease over the study period, but rather an alternation between positive and negative phases. After 2018, the common trend increased towards 2020, although it remained slightly negative in 2020 (−0.11).
Table 2 Key phases of the common primary production signal extracted by Dynamic Factor Analysis (DFA) for the five study areas in the central and southern eastern Adriatic Sea during the period 1998–2020. Positive DFA values indicate years with above-average common PP conditions, whereas negative values indicate years with below-average common PP conditions.
Tablica 2. Ključne faze zajedničkog signala primarne proizvodnje izdvojenog dinamičkom faktorskom analizom (DFA) za pet istraživanih područja u srednjem i južnom istočnom Jadranskom moru tijekom razdoblja 1998. – 2020. Pozitivne DFA vrijednosti označavaju godine s iznadprosječnim zajedničkim uvjetima PP-a, dok negativne vrijednosti označavaju godine s ispodprosječnim zajedničkim uvjetima PP-a.

Figure 4 Common interannual trend in satellite-derived primary production extracted by Dynamic Factor Analysis for the five study areas in the central and southern eastern Adriatic Sea during the period 1998–2020. Positive values indicate above-average common PP conditions, whereas negative values indicate below-average common PP conditions.
Slika 4. Zajednički međugodišnji trend primarne proizvodnje temeljene na satelitskim podacima izdvojen dinamičkom faktorskom analizom za pet istraživanih područja u srednjem i južnom istočnom Jadranu tijekom razdoblja 1998. – 2020. Pozitivne vrijednosti označavaju natprosječne zajedničke uvjete PP-a, dok negativne vrijednosti označavaju ispodprosječne zajedničke uvjete PP-a.
Source: Authors’ own work based on BICEP/NCEO data
Linear trend analysis was used to examine whether annual mean PP changed directionally at individual study areas during the 1998–2020 period (Table 3, Figure 5). Long-term trends were not uniform among study areas. A statistically significant negative trend was detected at Split Channel, with a trend coefficient of −0.99 mg C m⁻² d⁻¹ yr⁻¹ (R² = 0.20, p = 0.033). In contrast, Stončica showed a statistically significant positive trend, with a trend coefficient of +0.62 mg C m⁻² d⁻¹ yr⁻¹ (R² = 0.26, p = 0.013).
The remaining study areas did not show statistically significant long-term trends. Albania had a negative trend coefficient of −0.95 mg C m⁻² d⁻¹ yr⁻¹, but the trend was not significant (R² = 0.07, p = 0.238). Neretva showed the strongest negative coefficient, −1.39 mg C m⁻² d⁻¹ yr⁻¹, but this trend was also not significant (R² = 0.09, p = 0.165). Dubrovnik showed almost no directional change over the study period, with a trend coefficient of −0.09 mg C m⁻² d⁻¹ yr⁻¹ (R² = 0.00, p = 0.845).
Figure 5 illustrates these area-specific differences. Split Channel and Stončica showed opposite significant trends, whereas Albania, Neretva and Dubrovnik were characterized by interannual fluctuations without a statistically confirmed long-term direction. The relatively low R² values indicate that linear trends explained only a limited part of the total interannual variability, even at study areas where the trends were statistically significant.
Table 3 Linear trends in annual mean satellite-derived primary production at the five study areas during the period 1998–2020. Significant p-values are shown in bold and marked with an asterisk (*), with p < 0.05 considered statistically significant.
Tablica 3. Linearni trendovi prosječne godišnje primarne proizvodnje temeljene na satelitskim podacima na pet istraživanih područja tijekom razdoblja 1998. – 2020. Značajne p-vrijednosti prikazane su masnim slovima i označene zvjezdicom (*), pri čemu se p < 0,05 smatra statistički značajnim.
| Study area | Trend | R2 | p |
| Split Channel | -0.99 | 0.20 | 0.033* |
| Stončica | +0.62 | 0.26 | 0.013* |
| Neretva | -1.39 | 0.09 | 0.165 |
| Dubrovnik | -0.09 | 0.00 | 0.845 |
| Albania | -0.95 | 0.07 | 0.238 |

Figure 5 Long-term trends in annual mean satellite-derived primary production at the five study areas in the central and southern eastern Adriatic Sea during the period 1998–2020. Coloured lines show annual mean PP values, while dashed black lines indicate fitted linear trends. Values are expressed in mg C m⁻² d⁻¹.
Slika 5. Dugoročni trendovi prosječne godišnje primarne proizvodnje temeljene na satelitskim podacima na pet istraživanih područja u srednjem i južnom istočnom Jadranu tijekom razdoblja 1998. – 2020. Obojene linije prikazuju prosječne godišnje vrijednosti PP-a, dok isprekidane crne linije označavaju prilagođene linearne trendove. Vrijednosti su izražene u mg C m⁻² d⁻¹
Source: Authors’ own work based on BICEP/NCEO data
4. DISCUSSION / Rasprava
The results of this study indicate that satellite-derived PP variability across the analysed areas cannot be interpreted as a spatially uniform signal. This is evident from the significant differences among study areas and from the contrast between areas with higher mean PP, such as Albania, Split Channel and Neretva, and the lower-productivity study areas Dubrovnik and Stončica. The highest mean PP was recorded at Albania (450 mg C m⁻² d⁻¹), followed by Split Channel (422 mg C m⁻² d⁻¹) and Neretva (414 mg C m⁻² d⁻¹), whereas lower values were obtained for Dubrovnik (350 mg C m⁻² d⁻¹) and Stončica (328 mg C m⁻² d⁻¹). This spatial structure is consistent with the interpretation that PP variability may reflect a combination of local coastal forcing, regional Adriatic variability and broader Adriatic–Ionian/Mediterranean-scale processes, rather than a single controlling mechanism. Such an interpretation is consistent with the known heterogeneity of the Adriatic Sea, where productivity differs between shallow river-influenced coastal areas and more oligotrophic offshore waters [5].
A comparison with published PP values further supports this interpretation, although the agreement appears stronger for the general spatial pattern than for absolute values. At Stončica, the mean satellite-derived PP obtained in this study (328 mg C m⁻² d⁻¹) was higher than the corrected in situ regime values reported for this study area, including 128 mg C m⁻² d⁻¹ for 1997–2008, 251 mg C m⁻² d⁻¹ for 2008–2013 and 154 mg C m⁻² d⁻¹ after 2013 [15]. This suggests that satellite-derived PP may overestimate absolute productivity at this open-sea study area when compared with corrected in situ measurements. However, the obtained value still falls close to the broader historical range reported for Stončica and the middle Adriatic, where open-sea values are generally lower than coastal values [5]. Therefore, the satellite product seems to reproduce the expected relative spatial pattern, with higher PP in coastal and river-influenced areas and lower PP at the more open-sea study area, but exact quantitative agreement with in situ measurements should be interpreted with caution. Previous long-term studies in the middle Adriatic have shown that PP and related biological variability are characterized by substantial interannual and decadal changes rather than by a single persistent trend. Variability at Stončica and other long-term observation sites has been discussed in relation to changes in hydrographic conditions, atmospheric forcing, water-mass exchange and broader Adriatic–Ionian circulation patterns [11,14,15,17–22]. These studies provide an important context for the present results, particularly because the contrasting trends observed among the analysed areas suggest that long-term PP variability may reflect the combined influence of local coastal forcing and broader regional oceanographic processes.
At the local scale, the relatively high PP values at Albania, Split Channel and Neretva are consistent with stronger coastal and land-derived influences. This interpretation is particularly supported by study area Neretva, which had both relatively high mean PP and the highest coefficient of variation (31.5%), indicating that productivity in this area was not only elevated but also highly variable. Albania also showed the highest mean PP and the largest seasonal amplitude, suggesting a strong contrast between productive and low-productivity periods. In contrast, Stončica had the lowest mean PP and the smallest seasonal amplitude, which is consistent with its more open-sea character and weaker direct influence of terrestrial inputs. This pattern agrees with the broader understanding that coastal and semi-enclosed systems may respond strongly to nutrient inputs, hydrodynamics and residence time, while open-sea Adriatic areas are generally more oligotrophic [5,27].
The seasonal results further indicate that PP variability includes both a common regional component and area-specific differences. All analysed study areas showed the same basic seasonal structure, with a spring maximum and an August minimum. Maximum PP was recorded in May at Albania, Split Channel, Dubrovnik and Stončica, whereas Neretva reached its maximum earlier, in April. Minimum PP occurred in August at all study areas, with values ranging from 223 mg C m⁻² d⁻¹ at Stončica to 272 mg C m⁻² d⁻¹ at Albania. This common timing is consistent with a possible regional influence of seasonal stratification on PP across the analysed areas. The spring increase may reflect improved light conditions and nutrient availability after winter mixing, while the August minimum is consistent with stronger summer stratification and reduced nutrient supply to the euphotic layer. This interpretation agrees with the general functioning of oligotrophic and seasonally stratified Mediterranean systems [1,2,4,28]. The subsequent increase in PP in September may reflect the initial weakening of summer stratification, allowing greater nutrient supply to the euphotic layer while light availability remains relatively high. At coastal study areas, this seasonal response may additionally be influenced by freshwater and nutrient inputs. Because the September increase was observed across all analysed study areas, including the more offshore Stončica area, regional seasonal forcing may contribute to this pattern.
However, the magnitude of the seasonal cycle was not uniform among study areas. This is shown by the significant study area × season interaction and by the different seasonal amplitudes observed among study areas. Albania had the largest seasonal amplitude (407 mg C m⁻² d⁻¹), while Stončica had the smallest (264 mg C m⁻² d⁻¹). Neretva differed from the other study areas by reaching its maximum in April rather than May, possibly indicating an earlier response to riverine input.
The interannual results also indicate a combination of shared and area-specific variability. The strongest positive anomaly was recorded in 2010 at Albania (z = +3.20), and positive anomalies were also observed in the same year at Split Channel (z = +2.07) and Neretva (z = +1.64). This coincides with the strongest positive common signal extracted by DFA in 2010, when the DFA value reached 1.75. In contrast, reduced PP was particularly evident in 2018 at Neretva (z = −2.26) and Split Channel (z = −1.73), while the DFA also showed negative common conditions in 2017 and 2018. The response of the PP signal was not uniform. Stončica had its strongest positive anomaly in 2020 (z = +2.31), not in 2010, and its correlations with Neretva and Split Channel were weak. This indicates that the open-sea study area did not fully follow the same interannual pattern as the coastal and river-influenced study areas.
The broader Adriatic–Ionian context may nevertheless be relevant for interpreting part of the open-sea PP variability, particularly at Stončica. However, temperature, salinity, nutrient concentrations, freshwater input and circulation indices were not analysed directly in the present study. Consequently, the effects of BiOS, water-mass exchange, stratification and other hydroclimatic processes are discussed here as plausible oceanographic interpretations of the observed PP patterns rather than as statistically tested relationships. The BiOS mechanism modifies the pathways of Atlantic and Levantine waters, affects Adriatic thermohaline properties, and can influence nutrient redistribution and biogeochemical conditions in the Adriatic–Ionian system [20-22]. However, the PP anomalies identified in this study did not consistently coincide with the documented timing of Northern Ionian Gyre reversals. The strong positive common PP signal observed in 2010 preceded the circulation reversal reported around 2011, while no comparable PP response was evident around the 2006 reversal. This indicates that PP variability should not be interpreted as a direct response to individual circulation reversals. Instead, gradual changes in water-mass properties, nutrient redistribution, vertical mixing and stratification occurring during broader periods of Adriatic-Ionian circulation reorganization may influence PP before or around the timing of individual reversals. The observed temporal correspondence should therefore be regarded as a possible regional context for PP variability rather than evidence that circulation reversals directly controlled PP variability.
The significant positive trend at Stončica provides an additional indication that open-sea PP variability may be influenced by broader hydroclimatic processes rather than by direct local eutrophication. Stončica had the lowest mean PP and the smallest seasonal amplitude, but it was the only study area with a statistically significant positive trend. This combination is important: the study area is generally less productive and less seasonally variable, yet it shows a long-term increase. Such a pattern may be interpreted within the context of open-sea processes, including vertical mixing, nutrient redistribution, water-mass exchange, circulation-related variability and long-term changes in surface-water characteristics. However, these potential drivers were not tested directly in the present study. This interpretation is also consistent with previous work showing that long-term PP variability at Stončica includes regime-like shifts rather than a simple monotonic trend [15].
In contrast, Split Channel showed a different long-term pattern, with relatively high mean PP but a statistically significant negative trend. As an urban coastal area, this pattern may be associated with changes in local nutrient loading, wastewater management or other coastal environmental conditions; however, these potential drivers were not tested directly in the present study. Similarly, Albania and Neretva had relatively high mean PP, but their long-term trends were not statistically significant.
5. CONCLUSIONS / Zaključci
This study showed that satellite-derived PP across five selected areas of the central and southern eastern Adriatic Sea exhibited pronounced spatial, seasonal and interannual variability during the 1998–2020 period. The highest mean PP values were observed at Albania, Split Channel and Neretva, while lower values were recorded at Dubrovnik and Stončica. This spatial pattern indicates that coastal and river-influenced areas are generally more productive than the more open-sea study area, but also that productivity cannot be interpreted as a uniform regional signal.
A clear seasonal cycle was observed at all study areas, with a spring maximum and an August minimum. However, the magnitude of the seasonal cycle differed among study areas, as shown by the significant study area - season interaction. This pattern is consistent with a common seasonal component across the study areas, while the differences in seasonal amplitude indicate substantial area-specific variability. Neretva, with the highest coefficient of variation and an earlier seasonal maximum, showed a pattern consistent with stronger coastal and riverine influence.
Interannual variability revealed both shared and area-specific patterns. The DFA analysis identified a common regional signal, with the strongest positive phase in 2010 and negative phases in 2017–2018. However, individual study areas did not respond identically. Stončica showed a significant positive trend, whereas Split Channel showed a significant negative trend, while Albania, Neretva and Dubrovnik did not show statistically significant long-term trends. These results indicate that long-term PP changes differ among the analysed areas rather than representing a spatially uniform response within the study area.
Author contributions: Conceptualization, F.M. and G.J.M.; Methodology, Data Curation and Visualization, F.M. and N.G.; Formal Analysis and Writing – Original Draft Preparation, N.G.; Investigation, N.G., F.M., G.J.M. and M.K.; Writing – Review & Editing, N.G, G.J.M., F.M. and M.K.; Supervision, F.M.; Project Administration and Funding Acquisition, M.K.
Funding: This research was funded by the European Union (NextGenerationEU) under the Croatian Recovery and Resilience Plan 2021–2026 (NRRP), through the University of Split institutional project “Anthropogenic Impact on the Sustainability of Coastal Marine Ecosystems - AIME; IP-UNIST-52”, approved by the Ministry of Science, Education and Youth of the Republic of Croatia. Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor European Commission can be held responsible for them.
Conflict of interest: None.
Data Availability: The satellite-derived primary-production data used in this study are publicly available from the 9-km monthly primary-production dataset for 1998–2020 derived from ESA OC-CCI v4.2 products [24]. The geographical boundaries of the five extraction boxes are provided in Table 1, and the data were accessed through the OceanColour web portal (https://www.oceancolour.org). The R code used for data processing and statistical analyses is available from the corresponding author upon reasonable request.
Acknowledgments: The authors acknowledge the Ocean Colour Climate Change Initiative (OC-CCI) project for providing the satellite ocean-colour data used in this study.
Acknowledgement of AI or AI-assisted tools use: During the preparation of this work, the authors used ChatGPT for language editing, including grammar correction and clarity of expression. After using the tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the published article.
REFERENCEs / Literatura
[1] D’Ortenzio, F., & Ribera d’Alcalà, M. (2009). On the trophic regimes of the Mediterranean Sea: A satellite analysis. Biogeosciences, 6, 139-148.https://doi.org/10.5194/bg-6-139-2009
[2] Siokou-Frangou, I., Christaki, U., Mazzocchi, M. G., Montresor, M., Ribera d’Alcalà, M., Vaqué, D., & Zingone, A. (2010). Plankton in the open Mediterranean Sea: A review. Biogeosciences, 7, 1543-1586.https://doi.org/10.5194/bg-7-1543-2010
[3] Uitz, J., Stramski, D., Gentili, B., D’Ortenzio, F., & Claustre, H. (2012). Estimates of phytoplankton class-specific and total primary production in the Mediterranean Sea from satellite ocean color observations. Global Biogeochemical Cycles, 26, GB2024.https://doi.org/10.1029/2011GB004055
[4] Lazzari, P., Solidoro, C., Ibello, V., Salon, S., Teruzzi, A., Béranger, K., Colella, S., & Crise, A. (2012). Seasonal and inter-annual variability of plankton chlorophyll and primary production in the Mediterranean Sea: A modelling approach. Biogeosciences, 9, 217-233.https://doi.org/10.5194/bg-9-217-2012
[5] Matek, A., & Ljubešić, Z. (2024). Primary production research in the Adriatic Sea: A review. Acta Adriatica, 65, 205-221.https://doi.org/10.32582/aa.65.2.1
[6] Buljan, M., & Zore-Armanda, M. (1976). Oceanographic properties of the Adriatic Sea. Oceanography and Marine Biology: An Annual Review, 14, 11-98.
[7] Pucher-Petković, T. (1971). Recherches sur la production primaire et la densité des populations du phytoplancton en Adriatique moyenne (1962-1967) (Research on primary production and population density of phytoplankton in the Middle Adriatic (1962-1967)). In:21st CIESM congress, Monaco, 1968, 20 (3), Rapports et Procès-Verbaux des Réunions de la Commission Internationale pour l’Exploration Scientifique de la Mer Méditerranée, pp. 339-343.
[8] Pucher-Petković, T., & Marasović, I. (1988). Changes of productivity conditions in the open Middle Adriatic–eutrophication indicators. Pomorski zbornik, 26 (1), 585-593.
[9] Mifka, B., Telišman Prtenjak, M., Kuzmić, J., Čanković, M., Mateša, S., & Ciglenečki, I. (2022). Climatology of dust deposition in the Adriatic Sea; a possible impact on marine production. Journal of Geophysical Research: Atmospheres, 127 (7), e2021JD035783.https://doi.org/10.1029/2021JD035783
[10] Zore-Armanda, M. (1984). Hydrographic and productivity conditions of the Palagruža region in the Middle Adriatic. Acta Adriatica, 25 (1/2), 119-138.
[11] Matić, F., Grbec, B., & Morović, M. (2011). Indications of climate regime shifts in the Middle Adriatic Sea. Acta Adriatica, 52, 235-246.https://hrcak.srce.hr/148127
[12] Marasović, I., Viličić, D., & Ninčević, Ž. (1999). South Adriatic ecosystem: Interaction with the Mediterranean Sea. In P. Malanotte-Rizzoli & V. N. Eremeev (eds.), The Eastern Mediterranean as a laboratory basin for the assessment of contrasting ecosystems (pp. 383-405). Springer.https://doi.org/10.1007/978-94-011-4796-5_25
[13] Ninčević, Z., & Marasović, I. (1998). Chlorophyll and primary production of size fractionated phytoplankton in the Middle Adriatic Sea. In: 35th CIESM congress, Dubrovnik, 1998, 35 (2), Rapports et Procès-Verbaux des Réunions de la Commission Internationale pour l’Exploration Scientifique de la Mer Méditerranée, pp. 472-473.
[14] Marasović, I., Ninčević, Ž., Kušpilić, G., Marinović, S., & Marinov, S. (2005). Long-term changes of basic biological and chemical parameters at two stations in the middle Adriatic. Journal of Sea Research, 54 (1), 3-14.https://doi.org/10.1016/j.seares.2005.02.007
[15] Kovač, Ž., Platt, T., Ninčević Gladan, Ž., Morović, M., Sathyendranath, S., Raitsos, D. E., Grbec, B., Matić, F., & Veža, J. (2018). A 55-Year Time Series Station for Primary Production in the Adriatic Sea: Data Correction, Extraction of Photosynthesis Parameters and Regime Shifts. Remote Sensing, 10 (9), 1460.https://doi.org/10.3390/rs10091460
[16] Zore-Armanda, M., Stojanoski, L., & Vukadin, I. (1988). Time series of oceanographic parameters: Eutrophication of the open Adriatic waters. FAO Fisheries Report, 394, 71-77.
[17] Grbec, B., Morović, M., Beg Paklar, G., Kušpilić, G., Matijević, S., Matić, F., & Gladan, Ž. N. (2009). The relationship between the atmospheric variability and productivity in the Adriatic Sea area. Journal of the Marine Biological Association of the United Kingdom, 89 (8), 1549-1558.https://doi.org/10.1017/S0025315409000708
[18] Ninčević Gladan, Ž., Marasović, I., Grbec, B., Skejić, S., Bužančić, M., Kušpilić, G., Matijević, S., & Matić, F. (2010). Inter-decadal variability in phytoplankton community in the Middle Adriatic (Kaštela Bay) in relation to the North Atlantic Oscillation. Estuaries and Coasts, 33 (2), 376-383.https://doi.org/10.1007/s12237-009-9223-3
[19] Zore-Armanda, M. (1969). Water exchange between the Adriatic and the Eastern Mediterranean. Deep Sea Research, 16, 171-178.https://doi.org/10.1016/0011-7471(69)90072-2
[20] Civitarese, G., Gačić, M., Lipizer, M., & Eusebi Borzelli, G. L. (2010). On the impact of the Bimodal Oscillating System (BiOS) on the biogeochemistry and biology of the Adriatic and Ionian Seas (Eastern Mediterranean). Biogeosciences, 7, 3987-3997.https://doi.org/10.5194/bg-7-3987-2010
[21] Civitarese, G., Gačić, M., Batistić, M., Bensi, M., Cardin, V., Dulčić, J., Garić, R., & Menna, M. (2023). The BiOS mechanism: History, theory, implications. Progress in Oceanography, 216, 103056.https://doi.org/10.1016/j.pocean.2023.103056
[22] Eusebi Borzelli, G. L., & Carniel, S. (2023). A reconciling vision of the Adriatic-Ionian Bimodal Oscillating System. Scientific Reports, 13, 2334.https://doi.org/10.1038/s41598-023-29162-2
[23] Sathyendranath, S., Brewin, R. J. W., Brockmann, C., Brotas, V., Calton, B., Chuprin, A., Cipollini, P., Couto, A. B., Dingle, J., Doerffer, R., Donlon, C., Dowell, M., Farman, A., Grant, M., Groom, S., Horseman, A., Jackson, T., Krasemann, H., Lavender, S., ... & Platt, T. (2019). An ocean-colour time series for use in climate studies: The experience of the Ocean-Colour Climate Change Initiative (OC-CCI). Sensors, 19 (19), 4285.https://doi.org/10.3390/s19194285
[24] Kulk, G., Platt, T., Dingle, J., Jackson, T., Jönsson, B. F., Bouman, H. A., Babin, M., Brewin, R. J. W., Doblin, M., Estrada, M., Figueiras, F. G., Furuya, K., González-Benítez, N., Gudfinnsson, H. G., Gudmundsson, K., Huang, B., Isada, T., Kovač, Ž., Lutz, V. A., … & Sathyendranath, S. (2021). BICEP/NCEO: Monthly global marine phytoplankton primary production, between 1998–2020 at 9 km resolution (derived from the Ocean Colour Climate Change Initiative v4.2 dataset). NERC EDS Centre for Environmental Data Analysis. Retrieved September 24, 2026, fromhttps://doi.org/10.5285/69b2c9c6c4714517ba10dab3515e4ee6
[25] Zuur, A. F., Tuck, I. D., & Bailey, N. (2003). Dynamic factor analysis to estimate common trends in fisheries time series. Canadian Journal of Fisheries and Aquatic Sciences, 60 (5), 542-552.https://doi.org/10.1139/f03-030
[26] Holmes, E. E., Ward, E. J., & Wills, K. (2012). MARSS: Multivariate Autoregressive State-space Models for Analyzing Time-series Data. The R Journal, 4 (1), 11-19.https://doi.org/10.32614/RJ-2012-002
[27] Cloern, J. E., Foster, S. Q., & Kleckner, A. E. (2014). Phytoplankton primary production in the world's estuarine-coastal ecosystems, Biogeosciences, 11, 2477-2501.https://doi.org/10.5194/bg-11-2477-2014
[28] Moore, C. M., Mills, M. M., Arrigo, K. R., Berman-Frank, I., Bopp, L., Boyd, P. W., Galbraith, E. D., Geider, R. J., Guieu, C., Jaccard, S. L., Jickells, T. D., La Roche, J., Lenton, T. M., Mahowald, N. M., Marañón, E., Marinov, I., Moore, J. K., Nakatsuka, T., Oschlies, A., ... & Ulloa, O. (2013). Processes and patterns of oceanic nutrient limitation. Nature Geoscience, 6 (9), 701-710.https://doi.org/10.1038/ngeo1765
