Biophysical Consequences of a Relaxing Beaufort Gyre

A biophysical model shows that Beaufort Gyre (BG) intensification in 2004–2016 is followed by relaxation in 2017–2018, based on a BG variability index. BG intensification leads to enhanced downwelling in the central Canada Basin (CCB) and upwelling along the coast. In the CCB, enhanced downwelling reduces nutrients, thus lowering primary productivity (PP) and plankton biomass. Enhanced upwelling along the coast and in parts of the Chukchi shelf/slope increases nutrients, leading to elevated PP/biomass in the Pacific Arctic Ocean (PAO) outside of the CCB. The overall PAO PP/biomass is dominated by the shelf/slope response and thus increases during BG intensification. As the BG relaxes in 2017–2018, these processes largely reverse, with increasing PP/biomass in the CCB and decreasing PP/biomass in most of the shelf/slope regions. Because the shelf/slope regions are much more productive than the CCB, BG relaxation has the tendency to reduce the overall production in the PAO.


Introduction
Large-scale circulation in the Pacific Arctic Ocean (PAO; see Figure 1a for definition) is closely associated with the anticyclonic Beaufort Gyre (BG) located over the Canada Basin. Observations and model results indicate that the BG intensified in recent years; this intensification is associated with increased ocean velocity and freshwater content (FWC) in much of the PAO because of enhanced Ekman transport convergence and downwelling in the Canada Basin (CB) and enhanced upwelling in the Chukchi and Beaufort shelf and slope region (e.g., Giles et al., 2012;Krishfield et al., 2014;McPhee, 2013;Proshutinsky et al., 2009;Regan et al., 2019;Yang, 2009). However, there have been signs that the BG circulation began to stabilize in 2008 or 2009 (Armitage et al., 2017;Zhang et al., 2016).
The enhanced downwelling in the CB and upwelling in the Chukchi and Beaufort shelf and slope region likely have an impact on the planktonic ecosystem in the PAO. For example, based on observations taken in 2008, Coupel et al. (2015) reported that the CB during BG intensification showed a reduction in primary productivity (PP) because of a deepened nitracline (also see McLaughlin & Carmack, 2010), whereas areas with reduced freshening exhibited relatively high PP and phytoplankton biomass because of a shallower nitracline. While the impact of sea ice decline on Arctic PP and the planktonic ecosystem has been well recognized (e.g., Arrigo et al., 2008;Zhang et al., 2010;Jin et al., 2016), the impact of BG intensification and associated changes in the upwelling/downwelling pattern on PP and the planktonic ecosystem in the PAO has not been analyzed in a systematic manner on a decadal time scale. Even less known is the impact on the PAO planktonic ecosystem if significant BG relaxation occurs. Here a pan-Arctic biophysical model, the Biology-Ice-Ocean Modeling and Assimilation System (BIOMAS, Zhang et al., 2015), is used to examine changes in BG circulation and upper ocean physics during the period 1992-2018 and to assess how these changes affect the spatiotemporal variability of the planktonic ecosystem.

Brief Model Description
BIOMAS has been included in a number of community ecosystem model intercomparison studies (e.g., Jin et al., 2016;Lee et al., 2016). It consists of a sea ice model, an ocean circulation model, a pelagic biological model, and a sea ice algae model. The pelagic biological model has 11 components: two phytoplankton classes (diatoms and flagellates), three zooplankton classes (microzooplankton, copepods, and predatory zooplankton), dissolved and detrital particulate organic nitrogen, detrital particulate organic silica, nitrate, ammonium, and silicate (Zhang et al., 2015 algal components (diatoms and flagellates), with nitrate, ammonium, and silicate as limiting nutrients for ice algal growth in a 2-cm layer at the sea ice bottom. The ocean model is based on the Parallel Ocean Program (Smith et al., 1992). The sea ice model is adapted from the Pan-arctic Ice-Ocean Modeling and Assimilation System (Zhang & Rothrock, 2003), with melt ponds incorporated . BIOMAS assimilates satellite observations of sea ice concentration and sea surface temperature. It is forced by the National Oceanic and Atmospheric Administration 's Climate Forecast System (CFS) reanalysis data (Saha et al., 2010) over the period 1979-2018. More BIOMAS details are given in Zhang et al. (2015), Jin et al. (2016), and Lee et al. (2016), including model components, configuration, and initial conditions. Results over 1992-2018 are analyzed.

Model Evaluation and Results
In addition to assimilating satellite observations of sea ice concentration and sea surface temperature, BIOMAS sea ice velocity is calibrated with buoy drift data (http://iabp.apl.washington.edu) over the period 1992-2010, with a mean model ice speed bias of −6% and model-buoy speed correlation of 0.82. BIOMASsimulated isohaline (salinity = 31 psu) depth has a mean bias of~7 m (supporting information Figure S1) when compared to corresponding CTD (conductivity, temperature, and depth) observations in the central CB over the period 2003-2013 . The bias may be an indication that the model overestimates vertical mixing, in turn leading to a bias in PP. However, BIOMAS is able to capture most of the interannual variability of the CTD-derived isohaline depth, a measure of BG variability . When compared to available National Aeronautics and Space Administration IceBridge observations of sea ice thickness over the period 2012-2018 (supporting information Figure S2), BIOMAS shows a mean bias of 0.13 m (5%), with a model-observation correlation R = 0.64. A comparison of BIOMAS-simulated chlorophyll a (mg m − 2) integrated over the upper 100 m with corresponding observations compiled by the PacMARS project (Grebmeier et al., 2015(Grebmeier et al., ) over 1993(Grebmeier et al., -2012 indicates that most model results are relatively close to observations, especially in the region under study (supporting information Figure S3). However, the model also substantially overestimates or underestimates chlorophyll biomass in the upper 100 m at many observation locations, resulting in a large scatter in the plot of Figure S3. In particular, the model's inability to reproduce many of the observed large biomass values indicates the difficulty in comparing grid-cell-scale model results with point measurements. This was also reflected in the model-data intercomparison study of Lee et al. (2016) that used 21 biophysical models and a range of Arctic observations. This difficulty is often an important source of model bias, in addition to model uncertainties in forcing and parameterization. Overall, the model has a mean bias of −21.8 mg m −2 or 23% and the model-observation correlation is 0.35. A comparison of BIOMAS-simulated maximum chlorophyll a (mg m −3 ) in the water column with corresponding observations compiled by the PacMARS project over 1993-2012 (supporting information Figure S4) has similar features to those shown in Figure S3. ). However, the magnitudes of current speed and vorticity only decrease moderately in 2015, likely due to the long time scale of geostrophic current change in responding to changes in lateral density gradients (Johnson et al., 2018). In 2016, the average SLP remains low, so that the magnitude of vorticity continues to decrease in contrast to the current speed, which rebounds to some degree. Changes in SLP and the anticyclonic wind circulation affect upwelling and downwelling caused by the divergence (upwelling) or convergence (downwelling) of Ekman transport within the surface layer (Yang, 2009). The simulated Ekman upwelling velocity is calculated as w E = ∇×τ/ρf, where τ is ocean surface stress calculated following Martin et al. (2014), ρ is water density, and f is the Coriolis parameter. This method was shown by Zhong et al. (2018) to be highly correlated with similar methods that explicitly take into account the recent acceleration of geostrophic currents in the BG. The upwelling velocity is negative on average in the PAO (Figure 2b), indicating that downwelling dominates. Although the simulated w E is characterized by negative values (downwelling) over a large area of the CB, positive values (upwelling) occur over a  (Figure 1e), when compared to the 1992-2018 mean. In part of the Chukchi shelf and slope region, upwelling is increased. However, in 2017-2018, downwelling is reduced in most of the CB and upwelling is reduced along the coast (Figure 1f). In part of the Chukchi shelf and slope region, upwelling remains increased. Because the CB has a larger area than the coast, the magnitude of the PAO-averaged w E is reduced considerably (Figure 2b).
The model also shows decreasing sea ice extent, sea ice volume, and snow volume in the PAO during 1992-2018, in conjunction with increasing surface air temperature forcing (Figure 2d).  (Figure 2g) because of a deepened nitracline (e.g., Coupel et al., 2015;McLaughlin & Carmack, 2010). The simulated nitrate concentration is lower in the central CB and the East Siberian Sea and higher in the Chukchi Sea ( Figure 4a). As the BG intensifies during 2004-2016, nitrate concentration is reduced in the CB and increased in the Chukchi Sea and coastal areas in the Beaufort Sea (Figures 4a and 4b), owing to enhanced downwelling in the basin and upwelling along the coast and in part of the Chukchi shelf and slope region. Note that enhanced downwelling in the deep basin is associated with enhanced lateral fluxes (Ekman transport convergence) toward the basin from adjacent areas (shelf/slope regions). In the North Atlantic subtropical gyre, such enhanced lateral fluxes may lead to an increase in nutrients in the central area of the gyre (Doddridge & Marshall, 2018;Williams & Follows, 1998). In the BG, however, this increase does not happen because lateral nutrient fluxes are small due to the low level of simulated nitrate concentration in much the East Siberian Sea (Figure 4a) and high consumption of nutrients in the shelf/slope regions of the Chukchi and Beaufort seas as reflected in the elevated PP and plankton biomass (Figures 4e, 4h, and 4k).
When the BG relaxes in 2017-2018, the situation largely reverses, with the CB gaining nitrate and the Chukchi/Beaufort shelf region and the Beaufort slope region losing nitrate (Figure 4c). Meanwhile, nitrate in the Chukchi slope region is higher than the 1992-2018 mean, likely due to the fact that upwelling in 2017-2018 remains increased in part of the Chukchi shelf and slope region (Figure 1f Figure 2d). However, it is only weakly correlated with ocean temperature (0.34, p value = 0.08; Figure 2e). The insignificant correlation with ocean temperature indicates that the impact of changes in temperature on PP is limited in the PAO. The estimated temperature-dependent growth rate does not change significantly (<3%) in the temperature range of the region.
Increasing PP leads to generally increasing phytoplankton and zooplankton biomass before BG relaxation starts in 2017 (Figure 2h). Not surprisingly, the PP is well correlated with phytoplankton (R = 0.92) and zooplankton (0.86). During BG intensification in 2004-2016, PP and plankton biomass decrease in the central CB and increase elsewhere in the PAO, particularly in the shelf regions (Figures 4d,4e,4g,4h,4j,and 4k). This is likely because of the combined effect of the increase in light availability and changes in water temperature ( Figure 3e) and nitrate distribution (Figure 4b). During BG relaxation in 2017-2018, the model shows a large drop in average PP in 2017 and little rebound in 2018 (Figure 2g). The decreases in PP occur mostly in the areas away from the central CB, particularly in the shelf regions ( Figure 4f) mainly because of reduced nutrient availability (Figure 4c). In the central CB, PP tends to increase because of reduced downwelling. These changes in PP are translated to changes in plankton biomass, which show increases in the central CB and decreases elsewhere, including most of the shelf region (Figures 4i and 4l). These increases in plankton biomass in the central CB are even somewhat more pronounced than the increase in PP. The concurrent increase in both phytoplankton and zooplankton biomass suggests a lack of grazer control of phytoplankton biomass, which is supported by previous studies in this region (Campbell et al., 2009;Sherr et al., 2009).

Concluding Remarks
This model study examines the response of the planktonic ecosystem in the PAO to the BG intensification and relaxation, including regional (e.g., shelf vs. basin) responses. Consistent with various previous studies (e.g., Armitage et al., 2017;Giles et al., 2012;Krishfield et al., 2014;McPhee, 2013;Regan et al., 2019;Zhang et al., 2016;Zhong et al., 2018), the model shows strong BG intensification during roughly the period 2004-2016. The BG intensification results from more frequent occurrence of a strong Beaufort high-pressure atmospheric cell driving a strengthened anticyclonic wind and ocean circulation, as reflected in mostly positive values of the BGV index constructed based on the depth of a isohaline (S = 31 psu) in the central CB.
Model results further suggest that a significant relaxation of the BG has begun over the recent years 2017-2018. The BG relaxation is marked by negative values of the BGV index, with low SLP, ocean current speed and vorticity in the PAO.
The changes in BG circulation occurred at a time of sea ice decline, with a general increase in light availability at the ocean surface and in upper ocean temperature in the PAO, favorable for biological growth. The BG intensification with strengthened anticyclonic ocean circulation in 2004-2016 leads to enhanced Ekman transport convergence and downwelling in the central CB and enhanced upwelling along Alaskan coast.
In the central CB, the enhanced Ekman transport convergence of surface cold waters reduces water temperature and increases FWC, while the enhanced downwelling reduces nutrient availability because of a deepened nitracline (e.g., Coupel et al., 2015;McLaughlin & Carmack, 2010), resulting in reduced PP and plankton biomass. In most of the shelf and slope regions of the PAO, the enhanced Ekman transport convergence of surface waters toward the central CB and the general decreases in ice thickness cause an increase in water temperature. The enhanced upwelling along the coast and in some Chukchi shelf and slope regions causes an increase in nutrient availability, which, together with increasing light availability and water temperature, leads to elevated PP and plankton biomass in most of the PAO away from the central CB.
As the BG relaxes and the anticyclonic ocean circulation weakens in 2017-2018, the physical processes are largely reversed. The model simulates reduced Ekman transport convergence and downwelling in the central CB and reduced upwelling along the coast, when compared with 2004-2016. The reduced Ekman transport convergence causes a reduction in FWC in the central CB and an increase in most of other areas of the PAO. There is no significant drop in the total FWC in the upper 100 m of the PAO in 2017-2018 because of continued freshwater contribution from a thinning sea ice cover. In addition, the simulated water temperature continues to increase throughout the PAO, including the central CB. These physical changes in turn result in enhanced PP and plankton biomass in the central CB and reduced PP and biomass in most of other areas of the PAO including most of the shelf and slope regions, even though water temperature and PAR keep increasing there. Because the shelf and slope regions are much more productive than the central CB, BG relaxation in these recent years has the tendency to reduce the overall production in the PAO, as reflected in a drop in PP and biomass in 2017-2018 (Figures 2g and 2h). If the BG continues to relax in the future, we speculate that it may cause PP and plankton biomass to continue downward in the PAO. This could affect the CO 2 sink to the PAO and exacerbate warming of the PAO.
While model studies can shed light on changes in sea ice and upper ocean physics and how the changes may impact PP and the planktonic ecosystem in the PAO, it is essential to monitor the biophysical changes and possible interactions through satellite and in situ observations. Knowledge about the integrated system of sea ice, the upper ocean, and the planktonic ecosystem is still limited. It is necessary to learn more about the intertwining physical and biogeochemical processes in the PAO, for example, the close correlation between the BGV index and PP. This may be achieved by measuring changes in some of the key biophysical properties and linkages, such as changes in PP and the functioning of the planktonic ecosystem relative to changes in upwelling or downwelling and light or nutrient availability. These observations will enhance our understanding of the behavior of the integrated sea ice-ocean-biology system in a changing PAO with decreasing sea ice, increasing water temperature and FWC, and varying ocean circulation including BG intensification and relaxation. They will also help improve model representation of biophysical processes.