Seagoing data for the present study were collected between mid-2005 (Jul.) and mid-2015 at the outer Firth monitoring site ( Figure 1 ) over 40 voyages (predominately on NIWA’s RV Kaharoa but also using Western Work Boats Macy Gray and Star Keys). At approximately three-month intervals, seasonal voyages were undertaken to collect discrete biophysical data by deploying a Seabird SBE 911plus CTD (conductivity, temperature, depth) equipped with additional sensors for measuring Chl.a fluorescence (Seapoint Sensors Inc.), beam attenuation (CStar 660 nm, Wetlabs/Sea-Bird Scientific), and photosynthetically active radiation (PAR) using an LI-194 sensor (LI-COR Inc.). These measurements allowed for the profiling of water column properties. A rosette sampler with 12 x 10 L Niskin bottles fitted to the CTD captured water samples at six nominal depths (5, 10, 15, 20, 30 and 35 m) and processed according to methods detailed in Zeldis et al., 2022 and outlined below.

During each ship survey, a sub-surface mooring ( Supplementary Figure 1A ) equipped to measure water column biophysical properties was deployed, serviced, and subsequently re-deployed. Integrating Natural Fluorometers (INF-300, Biospherical Inc.) were positioned on the mooring to assess Chl.a (683 nm) and photosynthetically active radiation (PAR), within the upper surface mixed layer (~ 7 m) and deeper in the euphotic zone (~ 20 m) (see Chl.a methods below). The fluorescence and PAR sensors were kept free of biofouling using equipment that squirted saturated bromine solution on the sensor surfaces for 15 s every 3 h ( Supplementary Figures 1B, C ). Physical water column stratification was assessed using conductivity-temperature-depth-oxygen (CTDO) sensors (SBE 37 SMP-ODO, Sea-Bird Scientific) on the mooring within the mixed layer (~ 9 m) and near-bottom (~ 33 m) depths. All sensors sampled at 10-15 min intervals and daily means of data were used for comparison with seasonal ship survey and satellite observations.

Remotely sensed data obtained from NASA’s Moderate Resolution Imaging Spectrometer (MODIS), on the Earth Observing System (EOS) Aqua satellite (MODIS-Aqua), were utilized in this study (c.f. Figure 2 data-processing workflow outlined in Gall et al., 2022). To summarize, the raw data was received through direct broadcast at the NIWA Lauder satellite receiver in Otago. Calibration was performed, and ocean colour product algorithms applied to the data. Subsequently, the data were interpolated to a spatial resolution of 500 m on a transverse Mercator grid. For an extended analysis covering the period from mid-2002 to mid-2018, data were extracted from a 3 x 3, 9-pixel area centered around the study site. Medians of the nine-pixel area were calculated to optimize accuracy of the central pixel value which was within 5% in 99% of cases. This approach also effectively increased the number of observations throughout the timeseries by minimizing pixel failure caused by factors such as cloud cover, shadows, and breaking waves. As a result, there were approximately 20% more days available for observations.

CTD profiles were employed to assess various physical (temperature, salinity, density) and optical (PAR) properties and identify boundaries within the water column. These included: the surface mixed layer depth (z _m, m), where the density exceeded 0.125 kg m^-3 that of surface values (Levitus, 1982); the bottom mixed layer depth (z _bm, m₎, where the density reduced by 0.125 kg m^-3 that of bottom measures; the pycnocline width (y _pw, m), the distance (m) between layer depths (z _bm– z _m); the scalar light attenuation of PAR (K _oPAR), from the slope of the natural log transformed PAR profile; and the euphotic zone depth (z _eu), at which light is 1% of sub surface levels (Kirk, 2011).

In relatively clear waters satellite downwelling light attenuation of PAR (K _d PAR) can be estimated from a simple empirical relationship (Morel et al., 2007) to the K _d 490 band ratio algorithm product (Clark, 1997). In turbid waters a modified quasi analytical algorithm (QAA - Lee et al., 2005) from backscatter and the total absorption spectral coefficient is more accurate (Pinkerton, 2017). ‘Logistic-scaling’ of a backscatter product (b _bp 555, m^-1) was used in a blended procedure to estimate K _d PAR across Case 1 and Case 2 conditions (Pinkerton, 2017 and Gall et al., 2022). The NIWA long term monitoring site ( Figure 1 ) conformed to Case 1 conditions 98% of the time.

Statistical analysis and graphics were performed using R software (RStudio-Team, 2022), except for location maps (QGIS, 2022). Many of functions are include in base R, including linear model regression (lm) and statistical analysis functions from the Stats package (3.6.2). Standardized Major Axis (SMA) linear regression was used to minimize total least squares (TLS - orthogonal) rather than ordinary least squares (OLS) for bivariate comparisons using the ‘smatr 3’ package (Warton et al., 2012). This Model II regression method accounts for errors in both x and y variables to ensure statistically robust analysis. Figures were constructed using the grammar of graphics ‘ggplot2’ package (Wickham, 2011), customized for publication ready plots using the ‘ggpubr’ package (Kassambara, 2020). Missing daily data (mooring and satellite) were linearly interpolated between nearest values for seasonal statistics and analysis using the approx. function in the base R Stats package.

Seasonal in situ collections [mainly in early autumn (Mar.), mid-winter (Jul.), mid-spring (Oct.), and early summer (Dec.)], profiled a range of water column physical properties and strengths of vertical stratification ( Figure 2 ). Surface water column temperatures varied from mid-winter lows (about 13 °C) to highs (about 22 °C) in early autumn ( Figure 2A ). Surface salinities were generally lowest (<34.5 ppt) during mid-spring and highest (>35 ppt) during early autumn. The largest vertical difference in density stratification occurred in early summer ( Figure 2C ).

Consequently, surface mixed layer depth (z _m) and the pycnocline width (y _pw – in m not depth) were strongly and significantly seasonal ( Figures 3A, B and Table 1 ). Generally, mixed layer depth (z _m) was shallow between spring and summer (<10 m) and deep during mid-winter (>20 m), while pycnocline width (y _pw) was large (~ 20 m) during spring-summer and small (~ 5 m) during mid-winter ( Figure 3C ) - although at times exceptions occurred with multi-step water column density profiles evident ( Figure 2C ). The depth of the euphotic zone (z _eu) and its relationship to z _m and the depth to the base of the pycnocline (z _pb = z _m + y _pw) are relevant in terms of light supply for phytoplankton growth ( Figure 3D ). In most cases z _eu is deeper than z _m and comparable to z _pb, but at times deeper than both - particularly from early summer to early autumn. Seasonality in z _eu was not statistically significant and could attest to high interannual variability and low sampling numbers (n = 10, years). There is an indication of seasonality, with depths being generally shallower in winter and spring months, and deeper during the summer and autumn.

The time series of water column physical properties obtained from moored sensors and surface satellite products, collected on a quasi-daily basis ( Figure 4 ), provides significantly higher statistical power compared to the seasonal voyage snapshots taken every three months. These continuous measurements more reliably captured the seasonal periodicity on a monthly scale, as illustrated in Figure 5 . Satellite SST ( Figure 4A ) agreed well with mooring data ( Figure 6A , r ² = 0.91), although tended to have a slightly positive (higher) bias during summers and negative (lower) bias during winters, which may reflect different sampling depths (surface (<5 m) vs 7-10 m). Euphotic zone depth estimates from CTD and mooring data ( Figure 4C ) also agreed reasonably well ( Figures 6B–D , r ² = 0.55) although they differed at times. Satellite observations were like those from CTD profiles, with a small negative (lower) bias ( Figure 6C ), although regression and observation numbers are lower (r ² = 0.44, n = 17). Mooring z _eu contained a higher number of matchups with satellite data (n = 981), no bias, but lower agreement ( Figure 6D , r ² = 0.26), which is not surprising considering the different depth ranges of estimates (surface optical depth<5 m vs a moored sub-surface 7-20 m assessment). Complex water column stratification in physical properties and multi-stepped surface layers ( Figure 2 ) are likely to stratify light attenuating components across the water column such that at times, CTD, mooring and satellite observations sampled different surface layers.

All physical properties derived from mooring observations were seasonally significant in their periodicity regression models ( Figure 5 ). Mooring surface temperature was highest in late summer and lowest in late winter ( Figure 5A ). Surface salinity was highest in mid-autumn and lowest in mid to late spring ( Figure 5B ). Water column stratification was (generally) weakest during late autumn, and strongest during late spring ( Figure 5F ). There were obvious extremes and exceptions in patterns (e.g., lowest surface salinity water (<32 ppt) and largest salinity gradient (>1.5 ppt) during late winter 2008 ( Figure 4B ), coinciding with a flood event on 30-July-2018 ( Supplementary Figure 2B ). Other large flood events (e.g., Jan 2011 and 2012) also have low surface salinity, and both surface and deep salinities stayed generally lower during summers of these years compared with other years ( Figure 4B ). The lowest surface salinities observed during mid to late spring were typically delayed by approximately 1-2 months from the peak river flows that occurred in late winter ( Supplementary Figure 2B ). Euphotic zone depths in mooring records ranged from at least the maximum depth of the water column (about 40 m) to as shallow as 5-10 m in satellite records ( Figure 4C ). Strong z _eu seasonality was evident in mooring records across the year, with clearest waters (deepest euphotic zones) in mid-summer and more turbid waters (shallower euphotic zones) in mid-autumn; however, there was also considerable variation within months among years ( Figure 5C ). Shallower z _eu events (5-10 m) in satellite observations (2007, 2008 and 2009) than in mooring records ( Figure 4C ), coincided with near-surface salinity decreases ( Figure 4B ).

Seasonal CTD profiles of Chl.a (from fluorescence) spanning the period from 2005-2015, exhibited depth-dependent variations that coincide with physical stratification patterns ( Figure 2D ). During early summers, subsurface maxima layers (SCML) often formed at the base of the mixed layer or within the pycnocline, with high Chl.a values in the lower portion of the water column (HCL) compared to high Chl.a in the upper layer (HCU) - with an exception in 2006. Winter profiles generally had no SCML and associated with weak vertical physical stratification. Autumn and particularly spring, displayed a mix of higher Chl.a in upper and lower layers, with and without SCML at times.

Seasonal depth patterns are similar in the concentrations of phytoplankton biomass markers from water collections ( Figures 7 , 8 ). The distribution of Chl.a is like C _phyt above z _m but not below (white line, Figures 7A, C ). In general, C _phyt : Chl.a is higher above z _m (>50) than below (<50) ( Figure 7D ). High proportions of phaeopigment (Pha.a) highlight periods of senescence and decay ( Figure 7B ), which tended to occur in low Chl.a periods (i.e., generally opposite patterns to that of Chl.a - Figure 7A ). Further details on the size-partitioned of phytoplankton biomass are provided in Supplementary Figures 3–5 .

Applying averaging to water column sections above and below the depth of the mixed layer (z _m) provided a clearer depiction of seasonal patterns ( Figure 8 ). In early summer, above z _m, low Chl.a ( Figure 8A , ~0.6 mg m^-3) occurred with high C _phyt ( Figure 8C ), resulting in the highest C _phyt : Chl.a ratios ( Figure 8D , ~60). While the seasonality of the total biomass C _phyt : Chl.a ratio in surface waters was not significant (likely due to high interannual variability, and/or opposing size group effects), it was observed for pico-phytoplankton, with the highest ratios occurring in spring and summer ( Supplementary Figure 5H ). C _phyt was higher above z _m than below, especially from mid-winter to mid-spring ( Figure 8C ). Lowest concentrations of C _phyt occur below z _m in mid-winters (~40 mg m^-3), and above z _m in early summers (~60 mg m^-3). Chl.a and C _phyt seasonality patterns are similar below z _m. Seasonality was only significant for Chl.a above z _m, and Chl.a, Pha.a and C _phyt : Chl.a below z _m. The significant seasonal regression of the inverse ratio (Chl.a:C _phyt ), as presented in Table 1 , was utilized in the deep biomass CbPM_db-NPP modeling.

Phytoplankton biomass (Chl.a and C _phyt ) from mooring (7 and 20 m) and satellite surface (<5 m) data, ranged from about 0.3-7.5 mg Chl.a m^-3 and 50-500 mg C _phyt m^-3 ( Figure 9 ). Although satellite Chl.a was generally lower than mooring Chl.a ( Figure 9A ), seasonal patterns were similar, with high median concentrations during winter months and low concentrations during summer months ( Figures 10A, D ). Mooring Chl.a seasonality show a tendency for mid-autumn and late winter to early spring peaks, from a depression in early winter ( Figure 10A ), not seen in the more pronounced satellite seasonal cycle ( Figure 10D ). The ratio between upper (7 m) and lower (20 m) mooring Chl.a captured biomass-stratification seasonality, with high upper column Chl.a from autumn to spring (peaking in winter), to high lower column Chl.a from spring through summer (peaking in mid-summer) ( Figure 10B ; Table 1 ). Estimates of euphotic zone Chl.a from mooring surface Chl.a, and surface and deep Chl.a, were used to estimate a correction factor to apply to satellite surface Chl.a ( Figure 10C ; Table 1 ), to account for deep Chl.a in models.

Phytoplankton carbon (C _phyt ) from satellite estimates ( Figure 9C ) had similar average magnitude as near-surface C _phyt from seasonal surveys ( Figure 8C ) with finer temporal resolution. Satellite estimates had significant periodicity, peaking in mid-late winter, and low during summer ( Figure 10E ). Surface biomass C _phyt  : Chl.a ratios ( Figure 10F ) are similar in range and seasonality to those from seasonal surveys ( Figure 8D ), with significant (but low) seasonal periodicity. Biomass C _phyt : Chl.a ratio medians were generally low in early winter and high from spring through summer, although there was a large range of observations within months across the timeseries. Chl.a from mooring instruments agreed reasonably well with seasonal collections, with a small positive bias in near surface estimates ( Figure 11A , r ² = 0.27). Less bias occurred in deeper (20 m) estimates ( Figure 11B , r ² = 0.27). Chl.a from collections and from satellite (both from<5 m depth) had reasonable correlation (r² = 0.61) despite them having a limited number of satellite matchup days ( Figure 11C , n=26). Satellite and mooring Chl.a were significantly correlated but have low confidence (r ² = 0.18), with a small positive bias ( Figure 11D ). Correlations of C _phyt from seasonal collections with satellite estimates were weak (r² = 0.01, p=0.25), potentially due to the low number of observations ( Figure 11E , n = 16). For spring, summer, and winter there were suggestions of positive relationships, but not for autumn.

The Chl.a- based vertically generalized production model (VGPM) and the C _phyt based model (CbPM) supply two insights into assessing NPP from satellite surface data at the site ( Figures 12A, B ). Total annual NPP (integrated between yearly mid-winters (July) to span the production growth cycle) averaged between 220 g C m^-2 y^-1 and 161 g C m^-2 y^-1 for VGPM and CbPM respectively ( Table 2 ). Both models captured interannual variability of between about -10% to 10% of their overall annual climatological mean ( Table 2 ) and showed strong seasonality ( Figures 12C–F ). The VGPM annual mean results were about 36% higher than the CbPM, with less seasonality, mainly from much higher production rates in winter months and only slightly higher production rates during summer months. Both models assumed a degree of predictability between surface measure and areal z _eu estimates, which is clearly not the case at times at this site ( Figure 7 ; Supplementary Figures 7 , 8 ).

To account for deep biomass not detected in surface measures (and lack of predictability in areal euphotic zone estimates ( Supplementary Figures 7 , 8 ), we used the mooring Chl.a stratification (7m/20m) seasonality regression ( Table 1 ) to establish a seasonal correction factor (CF) for z _eu in the deep biomass VGPM_db model ( Figure 12A ), and a seasonal correction factor [Chl.a:C _phyt (z _bm)] ( Table 1 ) for deep carbon biomass in the CbPM_db model ( Figure 12B ). Accounting for deep biomass increased total annual NPP to between 231 g C m^-2 y^-1 (a 5% increase) and 191 g C m^-2 y^-1 (an 18% increase) for VGPM_db and CbPM_db respectively ( Table 2 ). Seasonally, accounting for deep biomass ranged from no change from mid-autumn to early spring, to largest increases during a late-autumn to an early summer peak for the VGPM (about 30%) and an early to mid-summer peak for the CbPM_db (about 33%).

Pinkerton et al., 2018 conducted an analysis of the accuracy of the MODIS-Aqua products utilized in our study region, offering valuable insights into satellite remote sensing of coastal water quality in the area. Furthermore, Gall et al., 2022 provides an overview of satellite remote sensing of coastal water quality in New Zealand, discussing both advances and challenges in this field, thus providing a broader context for our study.

It is worth noting that our moored instruments experienced tidal-driven depth variations of up to about 2.5 meters. The observed surface optical depths by the satellite, which were less than 5 meters, likely sampled shallower layers during periods of near-surface stratification compared to the moored instruments, which reached depths of around 7.5-10 meters. This phenomenon was evident at times in the CTD profiles depicted in Figure 2 and is supported by more detailed investigations into surface responses during discharge events within the Firth (Ø’Callaghan and Stevens, 2017).

Like other studies, our analysis reveals good agreement between sea surface temperature (SST) from satellite observations and in situ measurements. For instance, along the Korean Peninsula coast, Woo and Park, 2020 reported consistent results. Despite potential variations in assessment depths, we also found reasonable agreement in light attenuation (euphotic zone depth), with no evident seasonal bias. These findings align with comparisons made between satellite inherent optical properties (IOP) and field measurements in the China Sea, as demonstrated by Shang et al., 2011.

Our study revealed reasonable agreement among phytoplankton chlorophyll estimates obtained from 10 years of seasonal sampling voyages, moored instruments, and satellite observations, also without any notable seasonal bias ( Figure 11A ). However, there were discrepancies in the accuracy of the satellite and moored observations, with generally higher estimates from the mooring data, except for a few instances of low concentrations during summers. As highlighted above, differences may arise from variations in sampling depths, but also could be from methodological disparities between natural fluorescence measurements and satellite algorithms utilizing band ratio and absorption techniques. The mooring-based measurements of natural fluorescence were derived from daily averages which are known to include varying degrees of photochemical and non-photochemical quenching (Chamberlin et al., 1990). These quenching effects have the potential to influence the accuracy of the measurements.

Our satellite estimates of phytoplankton carbon, utilizing the empirical backscatter relationship of Graff et al., 2015 in the open ocean, demonstrated some agreement with in situ collections during mid-spring and early summer. However, in early autumn and mid-winter, the satellite estimates exhibited higher values. This higher bias coincided with periods when river stormflows, containing high loads of suspended mineral particles, were also seasonally high, potentially influenced the study site (see Supplementary Figure 1B ). It is well known that changes in the contribution of non-algal particles, including inorganic and organic particles (detritus), to the backscatter signal can introduce a positive bias (Stramski et al., 2004). Additionally, closer agreement and lower error has been observed when small pico-phytoplankton dominate biomass (Bellacicco et al., 2020), which supports our high pico-phytoplankton biomass and their ratios during mid-spring and early summer ( Supplementary Figures 8 , 9 ) and better agreement during these seasons than others.

Considering that the Case 1/Case 2 logistic scaling of backscatter for the site indicated Case 1 conditions 98% of the time, we conclude that our phytoplankton biomass products have similar accuracy to open ocean estimates most of the time. Further confidence in the estimates’ reasonableness is supported by similar biomass C _phyt : Chl.a ratios compared to in situ measurements (as shown in this study against others, Supplementary Figure 9 ). Additionally, our findings align with similar seasonal patterns observed in other satellite studies conducted in the South China Sea (Xu et al., 2021), surface coastal stations in the Black Sea (Stelmakh and Gorbunova, 2018, and Stelmakh, 2020), temperate Danish coastal waters (Jakobsen and Markager, 2016), and conservative agreement with data from the northwestern Atlantic (Sathyendranath et al., 2009).

The characterization of Chl.a stratification in coastal waters into different types, as presented by Zhao et al. (2019) for the German Bight, was observed at our study site. Our analysis of Chl.a biomass data obtained from mooring records yielded detailed monthly statistical distributions and boxplots that revealed consistent seasonal shifts in the distribution of Chl.a within the water column.

In early autumn to early spring, high Chl.a concentrations were predominantly found in the upper portion of the water column (HCU), while in late spring to late summer, they shifted to deeper regions either with a sub-surface maxima (SCML) or high Chl.a in the lower portion of the water column (HCL) to the bottom. These patterns were closely associated with the availability of dissolved macronutrients ( Supplementary Figure 6 ), primarily by inorganic nitrogen (DIN), including nitrate N and ammonium N, as reported by Zeldis et al., 2022 and Safi et al., 2022, modified by variations in light and physical stratification as discussed in Gall and Zeldis, 2011 and Zeldis et al., 2022. Surface DIN concentrations remained consistently low, below 1 μmol L^-1, throughout the seasons, except for early winter when they exceeded 1 μmol L^-1. This threshold indicates the potential for nitrogen limitation according to Eppley et al., 1969. Zeldis et al. (2022) and Safi et al. (2022) offer a comprehensive examination and insightful discussion of the alterations in DIN, dissolved reactive phosphorus (DRP), and silicate (DRSi) concentrations, shedding light on the nuanced connections between these changing ratios and shifts in phytoplankton taxonomic groups within the region. Despite the low DIN levels, the shallow stratified surface layers received sufficient light for biomass accumulation, peaking from late winter to early spring as light intensity increased. As DIN became depleted in the surface layer and temperature stratification strengthened and deepened from spring into summer, biomass accumulated in the pycnocline, progressing deeper into the lower water column (HCL) where DIN levels were not limiting (>1 μmol L^-1), and with sufficient light in clear surface waters (deep euphotic zones) and increasing light, would support biomass. Zhao et al., 2019 observed deep HCL layers during the decay phase of the spring bloom, associated with resuspension and erosion of the SCML, yet recognized the potential for photosynthetic activity below the pycnocline.

In a recent study spanning from March 2015 to October 2017 (c.f. Figure 9 ; Zeldis et al., 2022), the use of moored oxygen sensors at a depth of 35 m revealed a pronounced peak in gross primary productivity (GPP) in the deep layer during summers. This peak coincided with the highest positive values of apparent oxygen utilization, an indicator of net respiration, occurring prior to the breakdown of water column stratification in early autumn. As light levels decreased from mid-summer through autumn, deep biomass layers would become light limited. During the autumn to winter period, complete mixing of the water column occurred, as depicted in Figure 5 (and described by Zeldis et al., 2022 in their Figures 7 , 9 ), facilitating the replenishment of nutrient-depleted surface waters with macro-nutrients, completing the seasonal cycle.

The absence of cell senescence during mid-spring to early summer 3-month collections was supported by low phaeopigment ratios [Pha.a/(Chl.a + Pha.a)] that peaked in autumn and tended to be lowest in spring, albeit with interannual variability ( Figure 8B ). Safi et al., 2022 provided a detailed discussion on the development of a more heterotrophic system in early autumn, attributing the high ratios to the mineralization of spring primary production sometime during the mid-summer to early autumn period.

As well as nutrient supply and physical processes such as stratification, mixing and advective mechanisms influencing phytoplankton biomass and its stratification (Cullen et al., 2002), grazing, sedimentation, and cell lysis are also important mechanisms (Brussaard et al., 1995). Grazing pressure may also play a role in driving the formation and deepening of the deep Chl. a maxima (DCM) in the open ocean (Moeller et al., 2019), and previous studies in the inner Hauraki Gulf (Zeldis and Willis, 2015) showed increased rates of mesozooplankton grazing between spring and summer surveys that could add to biomass export to deeper layers.

An important distinction of deep Chl. a maxima DCM (SCML and/or HCL) is whether it is deep phytoplankton carbon (C _phyt ) biomass maxima (DBM) or deep photoacclimation maxima (DAM) - Cornec et al., 2021. Although we only documented SCML on some occasions in CTD profiles ( Figure 2D ), it is relevant to HCL at this deep estuary site (or conversely, shallow continental shelf site) in photoacclimation seasonality observed in our low C _phyt : Chl.a ratios in early summer ( Figure 8D ). As highlighted above, it was plausible that there was sufficient light (at times) to support deep HCL phytoplankton growth and maintenance, where increasing pigment concentration per cell in light limiting conditions is an adaption to harvesting more energy for photosynthesis (Geider et al., 1998 and MacIntyre et al., 2002). HCL at the Firth of Thames site coincided with clear surface waters and deep euphotic zones up to the mid-summer period of highest light intensities ( Figure 10B ).

In comparison to the study by Gall and Zeldis, 2011, our investigation provided a more comprehensive understanding of the seasonal dynamics of phytoplankton biomass in this deep estuary site. We detected significant seasonal variations in biomass ratios of pico-phytoplankton (mainly Prochlorococcus) in surface waters at our site, with minimum values in winter and maximum values in summer ( Supplementary Figures 6–7 ). Similar observations of larger sinusoidal variation were reported at the Hawaii Ocean Time-series (HOT) oligotrophic site in the North Pacific Gyre by Winn et al., 1995, as well as in satellite-based measurements of biomass ratios in the broader area around the HOT site (Westberry et al., 2016, c.f. Figure 2 ). Our research (Safi et al., 2022) also revealed an increase in micro-sized biomass during early summer and pronounced seasonality in deep biomass ratios, predominantly associated with micro-phytoplankton, originating from large centric diatoms and flagellates. The seasonal taxonomic shifts observed at our site may indicate that large motile phytoplankton with buoyancy control or flagella possess an advantage in exploiting nutrient gradients and light conditions during late spring and summer, as highlighted in studies by Sommer et al., 2017. It is well-known that diatoms, especially large centric diatoms, can rapidly adjust their sinking behaviors in response to changes in nutrient concentrations, as demonstrated by Du Clos et al., 2021.

The mean annual Chl.a based, and deep biomass corrected, VGPM-NPP (220-231 g C m^-2 y^-1 respectively) from July-2002 to July-2018 was higher than previous estimates (Gall and Zeldis, 2011) using a similar model from 1998 and 2000 (178 g C m^-2 y^-1). We attribute this difference to lower Chl.a estimates resulting from a larger averaging area that included deeper Hauraki Gulf waters ( Figure 1 vs c.f. Figure 1 ; Gall and Zeldis, 2011), rather than lower Chl.a at the Firth site.

In our study we used a fixed value of P _bopt (2.9 mg C (mg Chl.a)^-1 h^-1) determined in Gall and Zeldis (2011). It is important to note that P _bopt is known to be temperature dependent according to the findings of Harding et al. (2020). Although no seasonality was observed over the 13-23°C seasonal temperature range at the Firth monitoring site (Gall and Zeldis, 2011), this is likely due to the limited number of observations. If a greater availability of P _bopt data, like the findings in Harding et al. (2020) for Chesapeake Bay, had been present in our study, indicating a Pbopt range of 2.4-5.5 mg C (mg Chl.a)^-1 h^-1 across the 13-23 °C temperature range, it is likely that our analysis would have revealed more pronounced seasonality in VGPM-NPP.

The 25% lower CbPM-NPP models (161-191 g C m^-2 y^-1) are likely from lower growth estimates during the winter months, approximately half of those projected by the VGPM model. CbPM models are anticipated to better capture photoacclimation and physiological changes in C _phyt : Chl.a (Behrenfeld et al., 2005 and Westberry et al., 2008), rendering them more responsive to seasonal dynamics driven by factors such as nutrients and light, which have a significant influence on growth, supporting the larger range in seasonality of our CbPM models ( Figure 12 ).

Productivity models that incorporate depth-dependent responses of phytoplankton biomass and physiology to subsurface irradiance are generally considered more accurate in the deeper open ocean. However, these models only explain approximately 15% more variance in in situ NPP estimates compared to surface based VGPMs and CbPMs (Behrenfeld and Falkowski, 1997a and references therein). In our study, we incorporated stratified biomass and seasonality in deep Chl.a layers using mooring observations, which supported this finding, resulting in a modest increase of annual NPP estimates by 5% to 18% for VGPM and CbPM, respectively. Although one represents variance explained and the other represents an increase, we assume that these changes are influenced by similar factors. The largest increases in NPP occurred during the summer months (30-33%), coinciding with maximal deep biomass stratification. Moderate increases were observed in spring, while minimal or no increases were observed during autumn and winter.

Deep Chl.a maxima are commonly observed in various regions of the global ocean and as discussed above can arise due to either deep photoacclimation maxima or deep biomass maxima (Cornec et al., 2021), with the latter exerting a greater influence on net primary productivity (NPP). Previous investigations have highlighted the significant contribution of deep Chl.a maxima to water column NPP, particularly during specific seasons. For example, studies have shown that deep Chl.a maxima can account for up to 37% of NPP in the summer stratified conditions along the northern slope (30-80 m) of the Dogger Bank in the North Sea (Weston et al., 2005), and up to 60% during summers characterized by deep maxima in shallow coastal ecosystems (<50 m) of the northwestern Mediterranean (Guallar and Flos, 2017). Although these studies provide additional support for our findings, a quantitative assessment of in-water conditions is necessary to further validate these observations in our regions.

Different studies comparing bio-optical NPP models have shown varying preferences for a specific models in terms of their ability to accurately estimate NPP against in-water measurements (Saba et al., 2011; Silsbe et al., 2016; Regaudie-de-Gioux et al., 2019; Li et al., 2020). In our study, we employed the Chl.a based VGPM model, which has been previously used in the Hauraki Gulf continental shelf and Firth regions (Gall and Zeldis, 2011). This model is widely accepted due to its simplicity and reliable performance (Kahru et al., 2009). Additionally, we incorporated the carbon-growth-based CbPM model (Behrenfeld et al., 2005) as an alternative approach to evaluate the influence of biomass stratification in a straightforward manner.

Recent advancements have proposed mechanistic improvements for the carbon-growth-based CbPM model, utilizing a combined C _phyt , absorption, and fluorescence euphotic-resolved (CAFÉ) adaptive framework model (Silsbe et al., 2016). This enhanced model has demonstrated superior performance in capturing the highest variance and exhibiting the lowest bias when compared to other models on global in situ data. While the application of the CAFÉ model was not within the scope of this paper, it is planned for implementation as part of the NIWA-SCENZ (Sea’s Coasts and Estuaries NZ) framework for remote sensing of coastal water quality (Gall et al., 2022). Furthermore, it should be noted that accounting for deep biomass should also contribute to improving the accuracy of the CAFÉ model.

The net primary productivity (NPP) at the Firth monitoring site, adjusted for deep biomass, was estimated to be between 191-231 g C m^-2 y^-1, classifying it within the mesotrophic range (100-300 g C m^-2 y^-1) according to Nixon, 1995 coastal marine trophic classification. These productivity levels are similar to other estuarine-coastal ecosystems such as Delaware Bay and San Francisco Bay and close to the median of 181 g C m^-2 y^-1 in a worldwide compilation across 131 ecosystems in a comprehensive study by Cloern et al., 2014.

Our analysis of NPP data for the Firth site identified substantial interannual anomalies (-10% to 10%) ( Table 2 ). The monitoring site at the outer Firth should be equally impacted by nutrient inputs from both the seaward Hauraki Gulf to the north and its main river catchments to the south, as suggested in recent modelling scenarios (H. Macdonald, NIWA, pers. comm). This indicates that both sources play an important role in driving NPP variations at the site. In the Hauraki Gulf, the influence of wind-driven upwelling/downwelling events and mixing-driven nutrient fluxes along the shelf edge, has been described in studies by Sharples and Greig, 1998; Zeldis et al., 2004; Stevens et al., 2019). On the other hand, as highlighted previously, the nutrient concentrations within the Firth waters are predominantly influenced by catchment loading, as evidenced by the research of Zeldis and Swaney, 2018; Zeldis et al., 2022.

In their study, Safi et al. (2022) conducted trend analyses and utilized multivariate general additive modeling at our research site, revealing that micro-phytoplankton biomass demonstrated a primary response to elevated levels of dissolved inorganic nitrogen (DIN) and a secondary response to variations in physical stratification. These findings align with the observed seasonal responses in the biomass C _phyt : Chl.a ratio, which tend to be higher for smaller picophytoplankton when DIN concentrations are low. Similar patterns and correlations have been documented in more detailed investigations of nutrient conditions in the Baltic Sea transition zone (Jakobsen and Markager, 2016), where the ratios show a clear decrease with increasing nitrogen concentrations. The results of Safi et al. (2022) also suggest that physical stratification plays a role in driving variations in net primary productivity (NPP) at shorter temporal scales, such as seasonal fluctuations. Consequently, we anticipate that differences in the timing and intensity of nutrient loading and physical stratification across different years have contributed to the observed interannual variations in NPP at our monitoring site. However, further research is necessary to attain a comprehensive understanding and untangle the complexities associated with these processes.

Interannual variations in NPP, considering models adjusted for deep Chl.a (see Table 2 ), ranged from a low in 2008 (166-207 g C m^-2 y^-1) to a high in 2017 (215-259 g C m^-2 y^-1). Notably, the period between 2012 and 2016 exhibited consistently low annual NPP (~200 g C m^-2 y^-1), underscoring the challenges in interpreting long-term linear trends within specific multi-decade timeframes. An analysis of a longer time series of monthly mean Chl.a (used as a proxy for NPP) in the subtropical waters west of the North Island, spanning from 1997 to 2018 and combining data from SeaWiFS, MODIS-Aqua, and MERIS, revealed similar interannual variability. Their analysis highlighted non-linear multi-year fluctuations in Chl.a biomass and anomalies without any discernible trends (Pinkerton et al., 2019: c.f. Figures 3 , 4 , 5D ). Due to the small and varying trends observed over different time periods, as well as the need for a broader regional analysis encompassing the entire continental shelf from deep to shallow coastal margins, a comprehensive trend analysis for the region is being deferred to a future publication.

Chiswell et al., 2013 also conducted a study examining the climatology of surface Chl.a in the southwest Pacific Ocean surrounding New Zealand, focusing on the main water masses and fronts. Their research documented the phenology of phytoplankton bloom dynamics in the region. Specifically, they found that waters adjacent to the Hauraki Gulf typically exhibit short-lived mid-spring blooms, which occur when mixed layers shoal due to the relief of wind stress (Chiswell et al., 2013: c.f. Figure 5 ). These blooms occur for a brief period of less than one month and follow the replenishment of nutrients into surface waters during late autumn and winter mixing, providing sufficient light and nutrients to sustain phytoplankton growth. Gall and Zeldis, 2011 further explored these dynamics over the northeast continental shelf (at a water depth of 150 m) seaward of the Hauraki Gulf. They utilized SeaWIFS and moored surface INF Chl.a observations and detected a mid-late spring bloom in that area. In contrast, at the Firth monitoring site, where SeaWIFS and INF sampling were also conducted, a mixture of autumn, winter, and spring peaks was observed, but with earlier peaks occurring around winter and lower values during the summer period, consistent with the findings of the present study.

As the continental shelf nears the coastline, the depth available for biomass stratification decreases. We suggest that instead of the deeper chlorophyll maximum (DCM) or subsurface chlorophyll maximum (SCML) observed in open ocean waters, a high chlorophyll lower layer (HCL) forms in shallower waters, typically in less than 50 meters. This HCL is only found when surface waters are clear, when enough light can reach the lower layers for productivity. This transition is documented across a number of overseas studies in different regions, including: the Dogger Bank in the North Sea (c.f. Figure 2 , Fernand et al., 2013); the Australian continental shelf (c.f. Figure 8 , Mahjabin et al., 2020); the South Brazilian Bight (c.f. Figure 5 , Brandini et al., 2014); in NW Mediterranean coastal ecosystem (c.f. Figure 3 , Guallar and Flos, 2017); and summer glider transects across the continental shelf in the Mid-Atlantic Bight (Wright-Fairbanks et al., 2020). We have also observed similar transitions in Tasman and Golden Bays, New Zealand (Zeldis and Gall, 2005 and MacKenzie and Adamson, 2004), as well as across the northeastern New Zealand continental shelf into the Hauraki Gulf and Firth of Thames (c.f. Figure 2 ; Bury et al., 2012 and c.f. Figure 9 ; Safi et al., 2022). Additional investigation is necessary to understand, verify, and quantify the potential productivity of the high chlorophyll lower layer (HCL) and its correlation with areal net primary productivity (NPP).