The sampling site was situated in northern Marguerite Bay near Rothera Research Station and is part of the RaTS-program (Figure 1; 67.570°S 68.225°W). Year-round observations of biological and physical processes at this site (15 m depth) started in 1997 (Clarke et al., 2008; Venables et al., 2013; Rozema et al., 2017a). The RaTS site is situated in water ~520 m deep and ~2 km from the nearest shore. Our extensive summer campaign (Nov 2013–Mar 2014) was run in parallel to the RaTS program. CTD and fluorescence measurements to a maximum depth of ~95 m were obtained using a live-feed from the CTD package (SeaBird 19+) equipped with an in-line fluorometer (WS3S, WETLabs), turbidity sensor measuring at 700 nm (ECO NTU, WETLabs), and a spherical sensor (SPQA, LICOR) for Photosynthetically Active Radiation (PAR). Multiple fixed (2, 15, and 75 m) and one variable (depth of fluorescence maximum; F_max) depths were sampled for pigment analysis by HPLC, macronutrients (N, P, and Si) and maximum quantum yield of photosystem II (PSII). Phytoplankton production through ¹⁴Cincorporation was measured at F_max. Surface oxygen isotopes (0 m) used in this study were part of the RaTS program supplemented with a small number (n = 7) of samples collected at 2 m. Desired sampling frequency during summer was twice per week, and once per week for photosynthetic parameters with ¹⁴C incubations (P-E). The fluorescence profiles were calibrated against Chl-a measured by HPLC and collected at 15 m and F_max. All samples were collected by deploying a 12–24 L Niskin bottle with a hand or electric winch on a small boat. After collection, the samples were transported in dark and insulated boxes to the laboratories for further processing. Daily wind speed data from Rothera station was obtained from the British Antarctic Survey (http://www.antarctica.ac.uk/met/metlog/).

Year-round measurements of PAR were collected every half hour at Rothera station using a fixed sensor (SKP215, Quantum Skye, UK). These measurements were integrated to estimate the daily irradiance dose at the surface (PAR_daily, in μmol photons m⁻² s⁻¹). Attenuation in the water column (K_d in m⁻¹) was derived from linear regression of natural log-transformed, PAR-values as measured with the sensor on the CTD (Kirk, 1983). Estimates of K_d were based on the observations between 5 and 50 m unless PAR was below detection levels at depths shallower than 50 m. The first 5 m were excluded due to frequent shading by the sampling boat. Also, strong deviations in PAR within the water column during a cast, most likely due to changing cloud cover, resulted in a skewed estimate of K_d. Therefore, these casts were split and K_d was calculated by using the curve which had the majority of measurements. These estimates were linearly interpolated to obtain estimates for K_d on days between CTD casts and subsequently used to obtain daily light dose estimates for every depth for every half hour using (Equation 1):

where E (μmol quanta m⁻² s⁻¹) is the daily irradiance dose at depth z (m) (Kirk, 1983).

Mixed layer depth (MLD) was defined using the depth at which the Brunt–Väisälä frequency (BVF) was highest, as recommended recently for the WAP (Carvalho et al., 2017). Ocean Data View v4.7.8 was used to calculate the BVF (Schlitzer, 2016). Oxygen isotope samples (δ¹⁸O) were analyzed as described in Meredith et al. (2008). We used the δ¹⁸O and salinity data to derive quantitative estimates of the freshwater contributions from sea ice melt and meteoric water, which is possible since freshwater of meteoric origin is isotopically lighter than that originating from sea ice melt. The two sources of meteoric water are precipitation and glacial meltwater (Craig and Gordon, 1965). Simple budget calculations have shown that glacial discharge is the largest contributor to the meteoric water prevalence in the coastal WAP environment, however this glacial melt may be injected over a broader depth range than direct precipitation, which necessarily impacts the ocean directly at the surface (Meredith et al., 2017).

Nutrient samples were collected as described earlier (Bown et al., 2017) or as described below. A cross comparison of these two different collection methods during the 2012–2013 season did not yield major differences. Subsamples (2 × 5 mL) for macronutrient analysis were filtered (prefilter: 0.8 μm, membrane: 0.2 μm, Acrodisc Supor PF, Pall, US) and stored in the dark at −20°C (nitrate, nitrite, and phosphate) or 4°C (silicate). Analyses for samples collected using both collection methods were conducted using a Technicon TRAACS 800 Auto-analyzer at the Royal Netherlands Institute for Sea Research, The Netherlands (Bown et al., 2017). Nitrite and nitrate were summed and are presented as N_Tot.

Collection for and measurement of phytoplankton pigments were done using the standard RaTS protocol (Rozema et al., 2017a). In short, phytoplankton cells were collected on GF/F filters (Whatman, US) under low light and at ~2C using a mild vacuum (<0.2 mbar) generated with a water jet pump. After filtration, the filters were snap frozen in liquid nitrogen and stored at −80°C until analysis. For analysis, samples were freeze dried for 48 h prior to extraction in 90% acetone (v/v) at 4°C in the dark for another 48 h. The pigments were separated using a Waters 2,695 HPLC system equipped with a Zorbax Eclipse XDB-C8 column (3.5 μm particle size) as described by van Heukelem and Thomas (2001) and modified by Perl (2009). Pigments were manually identified using retention times and diode array spectroscopy (type 996, Waters, US) before quantification. Calibration of the system was performed using standards (DHI LAB PRODUCTS, Denmark) for chlorophyll c₃, peridinin, 19′-butanoyloxyfucoxanthin, fucoxanthin, neoxanthin, prasinoxanthin, 19′-hexanoyloxyfucoxanthin, alloxanthin, lutein, chlorophyll b, and Chl-a.

CHEMTAX (v1.95) was used to estimate the abundance of various phytoplankton groups (Mackey et al., 1996). This program employs a factor analysis and steepest descent algorithm to find the best fit using initial pigment ratios from previous studies for eight different phytoplankton classes, namely prasinophytes, chlorophytes, dinoflagellates, cryptophytes, two types of haptophytes, and two types of diatoms (Supplementary Table 1, Wright et al., 2009, 2010; Rozema et al., 2017a). Haptophytes were included as two separate groups since Phaeocystis, the most abundant haptophyte in our system, has shown great plasticity in cellular pigment ratios with iron and light availability (as reviewed in van Leeuwe et al., 2014). Also, two diatom groups were included to allow for differentiation in chlorophyll c₃ containing diatoms such as Proboscia spp. and Pseudonitzschia spp. and non-c₃ containing species (unpublished data).

A total of 139 samples were divided over 3 bins. We calculated from city-block distances and clustered the samples according to Ward's method as suggested by Latasa et al. (2010) using Past 2.17c (Hammer et al., 2001). This dendrogram suggested binning of the samples collected at 0–7, 8–49, and 50–75 m bins and the smallest bin size was 40. The initial ratios were as used previously with the RaTS data (Supplementary Table 1; Bown et al., 2017; Rozema et al., 2017a). Neoxanthin and prasinoxanthin were excluded for the 50–75 m bin. Further CHEMTAX settings were as described previously (Kozlowski et al., 2011; Rozema et al., 2017a). CHEMTAX was run 60 times using randomized pigment: Chl-a ratios (±35% of the initial matrix) to obtain the best possible result per bin (Supplementary Table 1; Wright et al., 2009). Both haptophyte and diatom groups as used with CHEMTAX were pooled and renamed to “pooled haptophytes” and “pooled diatoms,” respectively.

A quenching correction was applied to the upper water column fluorescence, as collected during the CTD casts, to allow for an accurate estimate of the water column biomass. This correction was based on linear regression of Chl-a concentration at 2 and 15 m or F_max (the sample collected at the shallowest depth was used) as estimated from HPLC analysis. If F_max was >>15 m (start summer campaign: late November 2013–Early December 2013) then Chl-a measured at F_max and 2 m were used. The slopes of these regressions were then used to estimate biomass in relation to the calibrated CTD fluorescence for the top layer. We present the corrected phytoplankton biomass estimates as Chl-a_corr. Additionally, standing phytoplankton stocks were calculated by integrating the first 80 m of the water column using Chl-a_corr. The depth of 80 m was chosen as this is generally in the layer of winter water, near the deepest sample for HPLC (75 m) and light does not penetrate this deep during summer (Clarke et al., 2008).

Triplicate measurements of maximal quantum yield of PSII (F_v/F_m) using a WATER-FT PAM (Walz, Germany) in dark-adapted subsamples (40 mL) allowed for the assessment of the photophysiological state of the phytoplankton community (reviewed in Maxwell and Johnson, 2000). F_v/F_m-values < 0.5 are generally indicative of exposure to stress due to e.g., exposure to excess PAR and/or ultraviolet radiation or nutrient limitation.

Dissolved Inorganic Carbon (DIC) concentrations were used in the calculation of primary production from ¹⁴C measurements. Samples for the analysis of DIC were collected on 8 sampling days (January 15, 25, 30, February 3, 10, 17, 25, and March 3, 2014) in borosilicate glass bottles and analyzed within 20 h of sampling as presented in Jones et al. (2017). DIC concentrations were determined by coulometric analysis according Johnson and Sieburth (1987) and Dickson et al. (2007) using a VINDTA 3C instrument (Versatile INstrument for the Determination of Total Alkalinity, Marianda, Germany). A power function was used to described the relationship between DIC and HPLC derived chlorophyll a concentrations (R² = 0.94, p < 0.001) to estimate DIC concentrations for the calculation of primary production on missing sampling days.

Water collected at the depth of F_max was used to estimate phytoplankton productivity. PP was estimated by constructing P-E curves (n = 16) based on a modified protocol from Lewis and Smith (1983). Seawater subsamples (20 mL) were incubated in a photosynthetron for 4 h under 21 different light intensities (8–713 μmol quanta m⁻² s⁻¹; Osram Powerball HCl-TT 250W in a Jolly symmetrical armature) at in situ (±0.5°C) temperature after the addition of 0.37 MBq ¹⁴C-bicarbonate (PerkinElmer, The Netherlands; >1.7 TBq mmol⁻¹) in clear 60 mL polyethylene terephthalate (PET) bottles (Qorpak, US). After incubation, the samples were filtered onto polycarbonate membrane filters (Millipore, US; 25 mm, 0.4 μm pore size) using a gentile vacuum (<0.1 bar) on a filtration carousel (Millipore, US). After filtration, the filters were acidified with fuming HCl for ~20 min to vent off unincorporated ¹⁴C-bicarbonate. Filters were placed into 20 mL PET scintillation vials (Wheaton, US) and supplemented with 10 mL scintillation cocktail (Ultima Gold, PerkinElmer, The Netherlands). For time zero activity, triplicate samples were prepared, filtered over polycarbonate membrane filters and immediately acidified. For total activity, 10 μL Ethanolamine was added to 100 μL ¹⁴C spiked subsamples in triplicate in 20 mL PET scintillation vials containing a clean polycarbonate membrane filter. Time zero activity and total activity were then treated identically to the other samples. All samples were vortexed and left for 2 days to ensure maximum dispersal of activity throughout the scintillation cocktail before counting activity on a Tri-Carb 2900TR (Packard Instrument Company, US).

The data were converted to DPM after correction for time zero activity (Steeman-Nielsen, 1952), and normalized to Chl-a, as derived from the HPLC analyses, and DIC (see above), and corrected for metabolic discrimination (5%) between ¹²C and ¹⁴C (Ærtebjerg-Nielsen and Bresta, 1984). These data were fitted by least-squares nonlinear regression using the model Equation 2 of Platt et al. (1980), modified by MacIntyre and Cullen (2005):

where P (μg C μg Chl-a⁻¹ h⁻¹) is the carbon fixation rate at irradiance E, P_s (μg C μg Chl-a⁻¹ h⁻¹) the light-saturated maximum of carbon fixation in the absence of photoinhibition, α [μg C μg Chl-a⁻¹ h⁻¹ (μmol quanta m⁻² s⁻¹)⁻¹] the initial slope of the P-E curve, β [μg C μg Chl-a⁻¹ h⁻¹ (μmol quanta m⁻² s⁻¹)⁻¹] is a measure for photoinhibition, and P₀ (μg C μg Chl-a⁻¹ h⁻¹) is the carbon uptake or release without any irradiance. These parameters were considered characteristic for the phytoplankton throughout the water column at the day of collection.

Two different methods of calculating phytoplankton carbon uptake were employed. The classical, mechanistic approach relies on the estimations of the P-E parameters, biomass estimates and total irradiance for days without P-E incubations (M_Platt; Platt et al., 1980). These estimates were obtained as discussed below and run through M_platt for every 1 m. Our second statistical approach method uses Random Forest (RF) regressions. RF is a machine learning method which constructs a large number of regression trees by randomly taking subsets of the data and input variables (Breiman, 2001). The major benefits of using such a tool over more traditional (mostly linear) models are: Regression trees allow for (1) non-linear relationships between input and output variables, (2) non-additive interactions between predictor and response variables and (3) correlations between different input variables (Breiman, 2001; Shi and Horvath, 2006). The relevant variables included were phytoplankton biomass, macronutrients, MLD, temperature, salinity, total irradiance and the relative abundance of the three most important phytoplankton groups, namely diatoms, haptophytes and cryptophytes (Behrenfeld and Falkowski, 1997; Arrigo et al., 1999; Dierssen and Smith, 2000; Vernet et al., 2008, 2012; Montes-Hugo et al., 2009; Piquet et al., 2011; Venables et al., 2013; Williams et al., 2016; Bown et al., 2017; Rozema et al., 2017a,b).

For the RF approach, we modeled the log of the PP estimates from the 16 P-E experiments (see above) using the depths for which all the aforementioned variables were available (n = 62). We constructed a model in which all the parameters (M_RFmax) were included and one reduced, optimized model (M_RFmin). During the construction we grew 2,000 trees, and four (M_RFmax) or two (M_RFmin) variables were used at every split. Variables least important in M_RFmax were removed individually until the fit of the new model no longer improved and this became M_RFmin. Partial dependence plots were generated for each variable per model to understand to the relationship with PP. These plots have the input parameter on the horizontal axis and give the range of impact of the input variable on PP. We discussed the importance of each variable to the models above but here we used the shape (and direction) of the curve to understand how this parameter influenced PP in our RF models. Partial dependence plots show the effect of the predictor on the response variable after accounting for the average effects of the other predictors. Therefore, it is the trend, not the actual values, that describes the dependence between both variables.

To assess the predictive power of both RF models, the data were split randomly into a training (80%) and test set (20%). This was repeated 10 times resulting in slightly different models of each type (M_RFmax and M_RFmin). These 10 slightly different models were subsequently tested against the randomly generated test sets using linear regression. We included the mean p and R²-values of these assessments. The “randomForest” package in R (v3.3.1) was used to build our models (Liaw and Wiener, 2002; R Core Team, 2016).

After the construction of the minimal RF model, we aimed to calculate the PP in the water column (0–80 m) over the full course of the season using M_platt and M_RFmin. For this, we needed to obtain full matrices for all parameters used in the models. We used anisotropic ordinary kriging to interpolate both in depth and time using the “gstat” package in R (Pebesma, 2004). First, we optimized the kriging function to best fit the autocorrelation per variable for the depth component using variograms. This function was then used for the interpolation through time. We chose a weighing factor of 0.3 for the time component to allow for the required sensibility needed during short, sudden events such as wind mixing. The matrix for PAR data was constructed as described above. Values for α, β, P₀, and P_s were linearly interpolated from days with P-E incubations. Thus, identical biomass and PAR estimates were used for all modeling approaches.

PP was calculated for every half hour, due to the resolution of the light measurements, and every depth (0–80 m) using M_platt and M_RFmin. The values were summed to obtain daily values. We chose to use a resolution of half an hour as light was highly variable throughout the day and our model included the potential for photoinhibition. Secondly, our RF model included the aforementioned parameters and was run using the constructed matrices. For light in M_RFmin used the summed daily irradiance for all depths. We summed all data for December, January and February as this is traditionally the most productive period and define this as the summer productivity estimates.

PAR availability in the water column was higher during spring (November) than during summer (Figure 2A). Also, PAR in the water column was lowest from December 15th to December 21th. This was related to high turbidity, increasing K_d to >0.30 m⁻¹ as well as low surface irradiance dose (December 23rd and 24th; Figures 2B,C). In contrast, surface irradiance dose peaked at the end of December before gradually decreasing. A brief and sharp increase of PAR at 15 m occurred on January 15thwhen K_d was low. Not until February 14th had PAR penetrated >50 μmol quanta m⁻² s⁻¹ beyond 15 m again and thereafter remained above this level (Figure 2A). The dynamics in K_d strongly coincided with standing phytoplankton biomass. The only major exception was after the short break (December 23rdand 24th) in our sampling effort due to hindrance of drifting sea ice sheets. During these casts, K_d was far higher (>0.30 m⁻¹) than observed at other times in our campaign.

Our field campaign started (November 26th) during a period where SSTs started to increase slightly while MLD and surface salinity decreased (Spring, Figures 3A,B, 4A). Above zero surface temperatures first occurred on and persisted onwards from December 6th(start of period 1). The MLD during the first half of December was shallow (1 m) while salinity kept decreasing. After a brief period of elevated wind speed and a deepening of the MLD (December 23rd), the water column surface salinity restabilized while the SST increased further (Figures 3A,B, 4A). Maximum SST at 1 m during this summer was 2.7°C on January 6th and 11th. The period of high SSTs was interrupted by an increase of wind speed and deepening of the MLD in mid-January ending period 1 and homogenizing water column temperature and salinity (Figure 4A). This increase in MLD was of a short duration as the water column restabilized within a week. Thereafter, surface temperature increased again with positive temperatures up to a depth of ~30 m and a stable surface salinity (~32.6). These conditions lasted until mid-February when strong winds occurred more frequently, salinity started to increase and SST decreased which resulted in the end of period 2. As in spring, the water column temperature dropped and stayed below 0°C and was more homogenous throughout the water column after February 10th, labeled as period 3.

Sea ice retreat from the preceding winter was late and the contribution of sea ice meltwater to the water column during summer was, although modest at all times, the highest (2.8%) recorded in the RaTS program (Figure 4B; Legge et al., 2017; Meredith et al., 2017). Oxygen isotope data from the surface and at 2 m depth showed an influx of sea ice melt after December 13th, continuing to at least the end of February. The contribution of glacial melt was higher than sea ice melt (mean 4.4%), but also more variable. After an initial increase, freshwater of meteoric origin peaked at 8.3% (also a record for the RaTS program) before dropping to just 3.1% on January 15th as an increase in wind speed mixed the water column (Figure 4A). Both freshwater components increased rapidly after the period of MLD deepening, after which the meteoric freshwater contribution at 2 m depth remained relatively stable at ~4% while the contribution of sea ice melt slowly decreased toward the end of summer.

Macronutrient concentrations at the end of spring were high with water column averages of 29.0, 1.8, and 83.5 μMfor N_Tot, phosphate, and silicate, respectively (Figure 5). These values were characteristic of the Antarctic winter water, as concentrations at 75 m were only slightly higher. During period 1, the vertical distribution of the macronutrients strongly followed the patterns observed in SSTs. Nutrient concentrations were generally lowest at the surface (2 m). Nutrient concentrations in the first ~15 m strongly decreased over a 10-day period after December 13th. On December 23, all macronutrients were virtually depleted with only 0.01 μM phosphate, 0.06 μM N_Tot, and 0.70 μM silicate left at 2 m depth. This period of depletion lasted until the end of period 1, when nutrients throughout the water column were mixed, restoring 2 m values to the levels similar to those observed before December 13th. Low nutrient concentrations were observed in period 2 with minimum values at 2 m depth on February 5th (13.6 μM silicate, 0.07 μM nitrate and no detectable nitriteor phosphate). This depletion of nutrients persisted to a depth of ~20 m. Period 3 was characterized by a gradual increase in surface nutrient levels.

Phytoplankton biomass accumulated during period 1 (Figure 3C). Relative abundances of diatoms increased at the cost of haptophytes during the start of this period (Figure 6). A sudden increase in phytoplankton biomass (up to 26 μg Chl-a l⁻¹), mainly consisting of small diatoms (Buma et al. unpublished microscopy data; Figure 6, Supplementary Figure 1), coincided with the occurrence of sheets of sea ice in the bay on December 14th. This biomass peak lasted until December 28th after which values dropped to ~5 μg Chl-a l⁻¹. Integrated phytoplankton stocks (0–80 m) followed a similar trend: a large increase between December 10th and 14th, and a steady decrease thereafter (Figure 7). Moreover, F_v/F_m in the top 10 m during the short, yet sharp increase in biomass was lower than expected for the in situ light intensities when compared to the rest of the summer (p < 0.001; Supplementary Figure 2). After the biomass peak, phytoplankton stocks remained stable until the end of period 1. During period 2, diatoms remained the most important taxonomic group although cryptophyte proportions increased (max. 40%) in the surface layer (2 m; Figures 4A, 6). Also, haptophytes increased at 15 m (~15%) until January 30th. Phytoplankton biomass was more evenly distributed over the first 30 m during this period, reaching a maximum concentration of 11 μg Chl-a l⁻¹ and suggested three short periods where biomass increased and decreased (Figure 7). The first was associated with the increase in haptophytes at 15 m, the second and third with diatoms. In period 3, relative abundances of haptophytes and cryptophytes were both elevated in comparison to the previous periods (max. 38 and 47%, respectively) and both phytoplankton groups co-occurred in the first 15 m (Table 2). Diatoms were more abundant at the deep fluorescence maxima. Phytoplankton biomass at 75 m was low and only increased during the two occasions where strong winds caused a deepening of the MLD. Relative group abundance matched those observed in the surface, albeit with a delay of ~30 days (Supplementary Figure 1).

Variability in photosynthetic parameters as obtained from the P-E experiments (n = 16) changed over the course of the summer. At the end of spring (November-early December), α, P_s, and P₀ were low in comparison to the summer values (Table 1). At the end of December, both α and P_s increased >50%. P_s appeared to reach an average of ~4 μg C μg Chl-a⁻¹ h⁻¹for the rest of the summer (January–March) while P₀ increased further to 0.68–0.85 C μg Chl-a⁻¹ h⁻¹. The α also increased further to 0.068–0.083 μg C μg Chl-a⁻¹ h⁻¹ [μmol quanta m⁻² s⁻¹]⁻¹. In contrast, β was very low and only increased to a maximum average of 0.004 μg C μg Chl-a⁻¹ h⁻¹ [μmol quanta m⁻² s⁻¹]⁻¹, suggesting that photoinhibition was of minor importance.

Using these P-E parameters, we calculated integrated PP for the water column on the days of the experiments (Figure 7). These values ranged from 92.5 (end of November) to 6,907.7 (end of December) mg C m⁻² d⁻¹. While the PP estimates from these experiments showed an increasing trend over the entire month of December, PP on December 23rd was substantially lower. The peak value of 6,907.7 mg C m⁻² d⁻¹occurred on December 30th after which PP decreased reaching a minimum of 438.2 mg C m⁻² d⁻¹on January 15th. PP during period 2 increased again to a maximum of 3,109.5 mg C m⁻² d⁻¹ and lasted until mid-February before stabilizing around 680 mg C m⁻² d⁻¹ during period 3.

To estimate the potential influence of phytoplankton species composition on PP we compared the log PP at discrete depths where we collected supplementary data (if PP > 0, n = 42) against the relative abundances of the three major phytoplankton groups (Table 2). It appeared that while PP and both diatoms and haptophytes approached a monotonic relation, cryptophytes and PP did not. Linear regressions for these three comparisons yielded r²-values of 0.01, 0.27, and 0.14 (p < 0.001) for cryptophytes, haptophytes and diatoms, respectively. Thus, we opted for spearman rank order analyses, which are less sensitive for non-linear trends (Table 2). We observed a negative relation between PP and haptophytes whereas a larger diatom fraction resulted in higher PP (p < 0.001). However, the relation between cryptophytes and PP was no longer significant (p = 0.36). The opposing trends between haptophytes and diatoms vs. PP were supported by correlations between the phytoplankton groups (Table 2). Haptophytes and cryptophytes were positively related to each other while both were negatively related to diatom abundance (p < 0.001).

Using M_RFmax with all the parameters included (p < 0.001, R² = 0.91), we quantified the contribution of the different input parameters on PP using a variable importance plot (Figure 8). Of these parameters, daily PAR was the most important, increasing the mean square error (MSE) by 50.4%. The other two large components were phytoplankton biomass (26.7%) and water temperature (21.6%). The three macronutrients (N_Tot, phosphate and silicate) showed a smaller importance (20.4–17.3%). Salinity (15.5%), MLD (13.6%) and the relative contributions of the major phytoplankton species (7.3–12.2%) influenced the MSE to a lesser extent. The contribution of some of the 12 included parameters was limited, possibly because some of the variability was already explained by a correlated and more important parameter. Therefore, we reduced our first model (M_RFmax) to obtain a more minimalistic model (M_RFmin) with only three parameters that explain 93% of the variability: daily irradiance dose throughout the water column, phytoplankton biomass and N_Tot (Figure 8). The second-best option explained only slightly less of the MSE and included seawater temperature instead of N_Tot. The inclusion of either seawater temperature or N_Tot resulted in highly similar minimal models. A linear correlation analysis on all the samples in which N_Tot was measured (excluding 75 m samples) revealed a strong relation between the two parameters (p < 0.001, r² = 0.67).

After the construction of our models, we analyzed how each of these parameters influenced PP at the depths where all auxiliary data (generally 2, 15, 75 m, and F_max) were available using partial dependence plots (Figure 9). The partial dependence plots for the two different models were highly similar and thus we present M_RFmax as it includes all variables, those for M_RFmin are included in the supporting information (Supplementary Figure 3). Light was the most important factor in both models, the partial dependence plot suggests that when light is not available, PP is strongly decreased. Also, when light is above a certain threshold (~3 mol quanta m⁻² d⁻¹), its effect on PP remains unchanged. Phytoplankton biomass showed a similar trend. PP did not increase further above ~7 μg Chl-a l⁻¹, although the support of this plateau is limited given that it is based on only 10–20% of the data. The three macronutrients had similar trends in relation to PP: the lower the concentration, the higher the PP until it reaches a plateau. Silicate and N_Tot dynamics suggest two distinct plateaux in its relation to PP, coinciding with the concentrations in the two different bloom periods (before and after January 15th). Seawater temperature showed a strong bimodal effect where 0°C was pivotal, likely related to the freeze of thaw cycle of glacial and sea ice. In contrast, a decreasing salinity and MLD coincided with a gradual increase in PP. Finally, the contribution of relative species abundances show a strong positive effect on the presence of diatoms as the effect on PP is strongly reduced below ~93% diatoms. Roughly 70% of our data have decreased diatoms abundances suggesting a strong support for this observation. As such, increasing fractions of cryptophytes and/or haptophytes are associated with decreased PP.

PP over the full course of the summer season was estimated using two different approaches, the traditional M_Platt and our best RF model, M_RFmin (Figures 7, 10). Dynamics of both PP estimates were very similar (Figure 10), and both methods produced near identical values for the period after January 15th (periods 2 and 3). While the patterns in PP for December and early January (period 1) were similar, the magnitude of PP differed greatly (Figure 10). At most, PP estimated in December was a factor 2.5 higher using M_Platt than when calculated using M_RFmin. On average, M_Platt was 40% higher during the first month of the summer reaching a maximum of 6,908 mg C m⁻² d⁻¹ on December 31st (Figure 7). PP in December peaked twice in both models, namely on December 15th and 31st. After January 15th, three steep increases in PP occurred on January 21st, 30th, and February 8–12th. The maxima of these periods were of the same magnitude at those calculated for December with the M_RFmin model. In contrast, using M_Platt PP in January and February was only half of the values observed in December. A second noteworthy difference between the two methods of estimation is to which depth PP occurs (Figure 10). However, these differences contribute little to the summed PP (Figure 10). Total summer PP (Dec–Feb) was estimated to be 214.4 g C m⁻² using M_Platt and 176.1 g C m⁻² using M_RFmin, or 2.68 and 1.96 g C m⁻² d⁻¹, respectively.

During period 1, the sea ice gradually melted and freshened the surface layer thereby reducing MLD and promoting diatom growth, as is typical for the traditional WAP ecosystem (Ducklow et al., 2012a, 2013). Moreover, a strong increase of turbidity indicated the marked presence of non-photosynthetic material, presumably microalgal aggregates released from the bottom of the sea ice (Figure 2B; Legge et al., 2017). Low F_v/F_m-values and relatively low values for α further suggests that the diatoms in these aggregates were low light acclimated which is typical for the sea ice environment (Table 1, Supplementary Figure 2, Arrigo et al., 2014). The difference in F_v/F_m between the late December community and the rest of the season suggests that the late December community was of a different origin, presumably the sea ice, confirming that period 1 could be deemed representative for the traditional coastal phytoplankton dynamics.

Classical features of the sea ice environment were the sole dominance of diatoms and depleted silicate stocks, the latter unexpected given the coastal nature of the site and never documented prior to the 2013–2104 summer (Clarke et al., 2008; Annett et al., 2017; Bown et al., 2017; Cassarino et al., 2017; Rozema et al., 2017b). Even though macronutrient concentrations in the water column during this period were low and uptake by phytoplankton was likely hampered, nutrients could still have been available through high bacterial remineralization, sloppy grazing, and/or viral lysis which peaks after the end of the bloom (Ducklow et al., 2012b; Brum et al., 2015). Additionally, the build-up of phytoplankton biomass could also have been ended by strong grazing, as is often observed at the sea ice edge, an attribute not covered in the RaTS program (Ross et al., 2008; Saba et al., 2014; Steinberg et al., 2015).

Brief wind mixing events marked the starts of both periods 2 and 3, two periods uninfluenced by sea ice cover. These mixing events reconfigured the water column, lower SST, and mixed the freshwater to a greater depth, thus changing the conditions for phytoplankton growth, as evident from changes in the community composition and PP. In period 2, the water column was restabilized by a fresh and shallow surface layer, which was populated by cryptophytes which are often associated with glacial meltwater (Moline et al., 2004; Saba et al., 2014). Below this layer, haptophytes and later large diatoms (mostly Proboscia spp.) thrived. Period 3 was characterized by increased wind mixing causing a more homogenous water column and deeper MLDs, typical for conditions at the end of summer and observed more frequently at the WAP (Clarke et al., 2008; Montes-Hugo et al., 2009; Venables and Meredith, 2014; Rozema et al., 2017b). Thereafter, the community consisted of large proportions of haptophytes, diatoms and cryptophytes (Figure 6).

First and foremost, our statistical model showed the importance of light availability and biomass in the promotion of PP. While completely expected, it does establish further confidence in our statistical model, which needed to establish these relations from the training set as opposed to the mechanistic model where these relations were already established.

Nutrient concentrations during summer can be highly variable and even undetectable near the surface (Figure 5; Clarke et al., 2008; Bown et al., 2017; Rozema et al., 2017b). Governing of PP (or biomass accumulation) by macronutrients in coastal Antarctic regions have been described previously and suggested a potential limitation of nitrogen (Alderkamp et al., 2012b; Henley et al., 2017; Rozema et al., 2017b). Moreover, the nutrient dynamics are strongly influenced by the source of meltwater. Glacial melt is low in macronutrients and organic matter and dilutes the surface layer (Henley et al., 2017). This opposed to sea ice melt, a prominent source of organic matter and bacteria to remineralized the nutrients (as reviewed in Fripiat et al., 2017). Additionally, deep (wind) mixing replenishes the nutrient concentrations of the euphotic zone. These three different components control MLD in the coastal WAP system and explains why nitrogen is a better predictor for PP than MLD (Meredith et al., 2008; Rozema et al., 2017b). The positive effects of a low salinity, above 0°C water temperature and shallow MLD on PP support the importance of glacial melt in shaping the PP. But in this study, these parameters did not contribute significantly in predicting PP (Figure 8).

Moreover, our partial dependence plots suggest negative relationships between nutrient concentrations and PP (Figure 9) when above 12, 0.7, and 42 μM for N_Tot, phosphate and silicate respectively and neutral when below these thresholds. We consider these negative relations an effect from the phytoplankton community becoming dominated by only one (or at least very few) species for which the conditions are optimal when the community is blooming (Annett et al., 2010). The blooming species, which contributes most to PP, takes up most of the macro nutrients and dominates. We suspect that the blooming species draws down nutrients to concentrations in line with their uptake kinetics and half-saturation constants and, when reaching these minimum concentrations, form the plateaus in the relations between macronutrients and PP (Huisman and Weissing, 1999; Timmermans et al., 2004).

When we consider all the parameters included in our M_RFmax model, the added benefit of including phytoplankton community composition is limited (Figure 8). However, this does not imply that Chl-a specific production does not vary between the groups in the natural assemblages. The random forest model selects the parameters that explain the statistically highest possible amount of the MSE of PP. Yet, some of the variability explained by the top parameters is also be explained by mechanisms related to lower ranking parameters. Thus, while adding those parameters would not improve the model this does not mean that they are not ecologically or mechanistically important. For example, literature suggests that haptophytes favor a more unstable water column and our results confirm such a relation (Arrigo et al., 1999; Rozema et al., 2017a). However, haptophyte relative abundance is suggested to be of only minor importance in relation to PP by our models, despite a large difference between photosynthesis strategies between haptophytes and e.g., diatoms (Kropuenske et al., 2009; Alderkamp et al., 2012a). However, as we know haptophytes are related to deeper MLDs, and therefore lower SSTs, higher macronutrient concentrations and more variable PAR (Kozlowski et al., 2011; Rozema et al., 2017a). Thus, most of the variability in PP explained by knowing the fraction of haptophytes is already captured by the aforementioned parameters that rank higher. As such, we conclude that adding species information would not significantly improve models of PP in the coastal system. Sufficient knowledge of the physical, (photo)chemical properties and total phytoplankton biomass in the water column will suffice for predicting PP.

Calculation of PP by both approaches (M_Platt and M_RFmin) was strongly similar except during the second half of December (Figure 10, period 1). The most important difference between these ~3 weeks and the rest of the season was the presence of melting sea ice. Estimates using M_Platt were higher (53%) during period 1 than those using M_RFmin. However, integrated biomass during period 1 was not exceptionally high when compared to period 2 or 3. In fact, it was lower than expected given the PP estimates from both models. Apparently, the high PP as measured using the P-E parameters did not result in an increase in total phytoplankton biomass. A difference in Chl-a specific production between taxonomic groups is unlikely to have been played a major role given the low contribution of taxonomic information to PP in M_RFmax (Figure 9). The large discrepancy between productivity and phytoplankton biomass accumulation could be due to a number of reasons. (1) Phytoplankton released from the sea ice are low light acclimated and were lowering their Chl-a:carbon ratio while acclimating to higher irradiance intensities in the water column, thus growth would remain undetected using Chl-a as we used Chl-a as a proxy (Sakshaug and Holm-Hansen, 1986). (2) Phytoplankton Chl-a:carbon ratios remain relatively unchanged, but the phytoplankton mainly replenish their micro- and macronutrient reserves before biomass accumulation by growth (Arrigo, 2005). These two explanations are supported by the lower F_v/F_m, possibly indicative of photo inhibition, and low α observed for this population (Table 1, Supplementary Figure 2). Additional explanations are related to ecosystems dynamics. They assume that the high productivity is translated into phytoplankton biomass, but there is also a large loss factor. (3) Sea ice occurrence is often associated with krill and other grazing which can efficiently lower phytoplankton biomass (Behrenfeld, 2010). (4) Viral lysis as after a period of biomass increase, viruses often lyse which include the bursting of phytoplankton cells and thus lowering biomass (Brum et al., 2015). Finally, the last and possibly most likely option (5) would be fast sedimentation of microalgal aggregates dominated by heavily silicified diatoms (Riebesell et al., 1991). The last three theories are supported by the observation of a very high turbidity in the surface layer (Figure 2B), not occurring elsewhere in the summer despite similar phytoplankton biomass concentrations. Our fieldwork cannot conclusively exclude any of the five hypotheses mentioned above or combinations thereof, that created the observed dynamics in phytoplankton biomass and proportions thereof to PP during period 1. Most likely it is a combination of multiple of the aforementioned explanations. Thus, while M_RFmin resembles mechanistic models that estimate PP in the marine environment, it would greatly benefit from the inclusions of parameters better explaining the phytoplankton dynamics at the sea ice edge. Also, the relation between N_Tot and PP does provide extra information for future efforts modeling PP in the coastal Antarctic ocean.

Our two estimates of 176.1 and 214.4 g C m⁻² per summer are in fair agreement with multiple different estimates of PP from previous years. Two recent estimates of net PP of 146 g C m⁻² yr⁻¹ for 2009–2010 and an estimate of 192 g C m⁻² yr⁻¹ for 2005–2006 for the RaTS site are comparable to ours (Weston et al., 2013; Henley et al., 2017). The latter estimate also employed P-E experiments with a very similar set up to ours. A slightly lower productivity of 99.4 g C m⁻² over a 5 month period was observed at a more northerly site near Palmer Station, generally considered less productive than northern Marguerite Bay (Stukel et al., 2015; Rozema et al., 2017a). Furthermore, our estimates are in broad agreement with a multi-year (1995–2006) analysis of PP using ¹⁴C along the WAP (Vernet et al., 2008). Another study based on the 2007–2008 summer but also covering the large WAP region estimate a lower productivity (78 g C m⁻² over 4 months) although this estimate is from a low biomass summer (Huang et al., 2012; Rozema et al., 2017a). Another estimate derived from calculating NCP from silica isotopes conducted at the RaTS, has calculated a maximum of only 56.4 g C m⁻² over 4 months in 2009–2010 (Annett et al., 2017). Such a large discrepancy might be related to the bias of a silica based method to diatoms as haptophytes and cryptophytes also contribute to the carbon drawdown (Figure 9) and diatom abundance was low in the 2009–2010 seasons (Rozema et al., 2017a).

Our data underline the importance of sea ice for phytoplankton biomass accumulation and PP in the coastal seas. Thus, a decrease in sea ice, as has been observed for the WAP since the start of the satellite record, will have substantial impact on the coastal marine system (period 1; Stammerjohn et al., 2008). Also, the diatoms associated with the sea ice and the related meltwater layer are pivotal in the food web due to the importance in the krill life cycle (Flores et al., 2012). While the loss of sea ice would result in a decrease in phytoplankton, the influence of meltwater from glaciers is rising (period 2; Moline et al., 2004; Montes-Hugo et al., 2009). While this meteoric meltwater also promotes stratification, which in turn promotes PP, the timing of the establishment of glacial meltwater layers and the phytoplankton communities within are different potentially affecting the entire food web (Moline et al., 2004; Montes-Hugo et al., 2009; Mendes et al., 2013; Saba et al., 2014; Rozema et al., 2017a). While most PP models are designed for ocean regions, the MLD at the coastal ocean is established differently than in the open ocean due to the proximity of glaciers thus are generally not adequately included in these models.

Our minimal model also underlines that dissolved nitrogen and seawater temperature can be used to estimate PP. While we do not infer a causal relation between PP and these parameters, we observe strong relation to PP. Our RF does not fully explain the dynamics in PP, namely only 93%, and we consider this is related to the dynamics occurring in December. The strong reduction in phytoplankton biomass despite a high PP suggests another source which influences phytoplankton stocks, as discussed above. We consider the counterintuitive relation between N_Tot and PP evidence of important underlying processes not captured by our accurate yet simple model. Thus, a future improvement would be the measurement and incorporation of these undescribed processes into a model. Also, the near binary effect of above and below 0°C sea surface temperatures only indicates the importance of meltwater and mixing to the winter water layer. No relation to PP was observed with increasing seawater temperature while the strength of the relation between PP and salinity declines strongly when approaching lower values (<32.5).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer DCLS and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.