Development of adjoint-based ocean state estimation for the Amundsen and Bellingshausen Seas and ice shelf cavities using MITgcm/ECCO (66j)

. The Antarctic coastal ocean impacts sea level rise, deep-ocean circulation, marine ecosystems, and the global carbon cycle. To better describe and understand these processes and their variability, it is necessary to combine the sparse available observations with best-possible numerical descriptions of ocean circulation. In particular, high ice-shelf melting rates in the Amundsen Sea have attracted many observational campaigns and we now have some limited oceanographic data that capture seasonal and interannual variability during the past decade. One method to combine observations with numerical models that 5 can maximize the information extracted from the sparse observations is the adjoint method, aka 4D-Var, as developed and implemented for global ocean state estimation by the Estimating the Circulation and Climate of the Ocean (ECCO) project. Here, for the ﬁrst time, we apply the adjoint-model estimation method to a regional conﬁguration of the Amundsen and Bellingshausen Seas, Antarctica, including explicit representation of sub-ice shelf cavities. We utilize observations available during 2010–2014, including ship-based and seal-tagged CTD measurements, moorings, and satellite sea-ice concentration estimates. 10 After 20 iterations of the adjoint-method minimization algorithm, the cost function, here deﬁned as a sum of weighted model-data difference, is reduced by 65% relative to the baseline simulation by adjusting initial conditions, atmospheric forcing, and vertical diffusivity. The sea-ice and ocean components of the cost function are reduced by 59% and 70%, respectively. characteristics and (2) intrusions of modiﬁed Circumpolar Deep Water (mCDW) towards the Pine Island Glacier. Sensitivity experiments show that ∼ 40% and ∼ 10% of 15 improvements in sea ice and ocean state, respectively, can be attributed to the adjustment of air temperature and wind. This study is a preliminary demonstration of adjoint-method optimization with explicit representation of ice-shelf cavity circulation. Despite the 65% cost reduction, substantial model-data discrepancies remain, in particular with annual and interannual variability observed by moorings in front of the Pine Island Ice Shelf. We list a series of possible causes for these residuals, including limitations of the model, the optimization methodology, and observational sampling. In particular, we hypothesize that 20 residuals could be further reduced if the model could more accurately represent sea-ice concentration and coastal polynyas.


1
The ice shelves and glaciers in the Amundsen Sea (AS) and Bellingshausen Sea (BS) are melting and thinning rapidly with consequences for global sea level rise and changes in ocean circulation and the global carbon cycle (e.g., Arrigo et al., 2008;Pritchard et al., 2012;Paolo et al., 2015;Bronselaer et al., 2018;Rignot et al., 2019). Basal melting of these ice shelves is 25 caused by warm modified Circumpolar Deep Water (mCDW, 0.5-1.5 • C), which intrudes onto the continental shelf toward the ice shelf cavities following submarine glacial troughs ( Fig. 1) (e.g., Jacobs et al., 1996;Walker et al., 2007;Jacobs et al., 2011;Nakayama et al., 2013;Walker et al., 2013;Dutrieux et al., 2014). For this reason, multiple oceanographic observational campaigns have been collected by the international community to understand the mechanism of mCDW intrusions onto the AS continental shelf and towards ice shelf cavities. As part of these efforts, we now have some limited oceanographic data that 30 capture seasonal and interannual variability during the past decade (e.g., Jacobs et al., 2011;Nakayama et al., 2013;Dutrieux et al., 2014;Heywood et al., 2016;Kim et al., 2017;Webber et al., 2017;Mallett et al., 2018).
Recent observations as well as modeling studies reveal that mCDW pathways, ice shelf-ocean interaction, the thermocline depth, and ocean bathymetry below Pine Island Ice Shelf (PIIS) are important for controlling the PIIS melt rate (e.g., Schodlok et al., 2012;Dutrieux et al., 2014;De Rydt et al., 2014;St-Laurent et al., 2015;Dinniman et al., 2016;Jourdain et al., 2017; 35 Kimura et al., 2017;Webber et al., 2019). The thermocline depth was ∼200 m deeper in 2012 compared to other years (e.g., 1994, 2007, 2009, see Fig. 2A in Dutrieux et al. (2014), which reduced the PIIS melt by ∼ 50%. After 2012, the thermocline shoaled by 200m returning to its more commonly observed depth of ∼350 m (Webber et al., 2017). It is suggested that this thermocline variability was caused by changes in local and remote surface wind and buoyancy forcing (Dutrieux et al., 2014;Webber et al., 2017). 40 To better describe and understand these processes and their variability, it is necessary to combine the sparse available observations with best-possible numerical representations of ocean circulation. One method to combine observations with numerical models that can maximize the information extracted from the sparse observations is the adjoint method, also known as 4-Dimensional Variational assimilation (4D-Var), as developed and implemented for global ocean state estimation by the Estimating the Circulation and Climate of the Ocean (ECCO) project. To date, the ECCO project has produced ocean state es-45 timates based on Circum-Antarctic or global model configurations (e.g., Mazloff et al., 2010;Forget et al., 2015;Zhang et al., 2018;Fukumori et al., 2020). Employing the adjoint model produced by automatic differentiation (Giering and Kaminski, 1998), aka algorithmic differentiation, and utilizing temporally-varying oceanographic observations, these ocean state estimates are capable of simulating the large-scale evolution of the Southern Ocean consistent with the available observations. Many observational and modeling studies have been conducted to understand Southern Ocean gyre dynamics, subsurface 50 ocean circulation, the southern shift of various fronts around Antarctica, etc. (e.g., Gille et al., 2016;Jones et al., 2016;Tamsitt et al., 2017;Nakayama et al., 2018;Roach and Speer, 2019;Jones et al., 2020). However, despite the importance of Antarctic coastal regions for global climate, existing models fail to accurately reproduce the sparse available observations, likely owing to the difficulty in simulating Antarctic continental shelf regions and sub-ice-shelf-cavity processes (Mazloff et al., 2010;Tim-Kusahara, 2020).
For other regions of the globe, ocean state estimates based on regional configurations have been successfully developed during the past decades, achieving good model-data agreement and leading to understanding of reginal processes (Fenty and Heimbach, 2013b,a;Verdy et al., 2014;Rudnick et al., 2015;Nguyen et al., 2020;Verdy et al., 2017;Vinogradova et al., 2014). For Antarctic coastal regions, however, the only previous attempt to constrain a model with observations was the study 60 of Nakayama et al. (2017), which used a low-dimensional estimation approach based on the computation of model Green's functions. Here we aim to extend the study of Nakayama et al. (2017) by employing the adjoint method, which permits a larger number of higher dimension control variables than the Green's functions approach. The objective is to obtain a closer fit to the available observations than what was achieved in Nakayama et al. (2017). This objective is challenging due to several difficulties including (1) polar specific processes (ice shelf and sea ice) are highly nonlinear and (2) observational data is 65 limited. The groundwork for making adjoint-method optimization possible in the presence of ice shelf cavities was laid out in the study of Heimbach and Losch (2012), who obtained adjoint sensitivities of sub-ice shelf melt rates to ocean circulation under Pine Island Ice Shelf, West Antarctica.
In this study, we present our attempt at the development of Amundsen-Bellingshausen Seas ocean state estimates by employing the adjoint-model-based data assimilation method developed by ECCO for regional and global ocean state estimation 70 (Mazloff et al., 2010;Forget et al., 2015;Zhang et al., 2018;Fukumori et al., 2020). We focus on the years 2010-2014 when oceanographic observations were collected frequently and the largest interannual variability has been observed (Dutrieux et al., 2014;Webber et al., 2017). Our simulations are carried out for a subregion of the global 1/3 • ECCO solution, aka ECCO LLC270 . Using the ECCO LLC270 solution both provides lateral boundary conditions for this study as well as enabling this work to be a stepping stone towards improved representation of ice-shelf-ocean interactions in ECCO 75 global-ocean retrospective analyses. We note, however, that the LLC270 horizontal and vertical resolutions are insufficient to resolve critical ocean and ice-shelf processes, e.g., eddy transport and mean-flow-topography interactions. Hence, these subgrid-scale processes need to be parameterized and adjusted.
2 Data and methods

80
In the Amundsen Sea, oceanographic observational campaigns were carried out in (Nakayama et al., 2013Dutrieux et al., 2014;Heywood et al., 2016;Kim et al., 2017). Several mooring observations were also obtained, with the moorings at the PIIS front capturing the largest interannual variability observed in the region between 2009(Dutrieux et al., 2014Webber et al., 2017). We also utilize seal-tagged CTD observations obtained in 2014, which contain over 10,000 profiles between February and November (Heywood et al., 2016). In the central part of the Bellingshausen Sea, no oceano-85 graphic observations were collected between 2010-2014. Recently, we have become aware that seal-tagged CTD observations, mostly in the Bellingshausen Sea (Roquet et al., 2013;Zhang et al., 2016), are available for inclusion in future studies. For the Antarctic peninsula region, oceanographic observations were collected by the Palmer Antarctic Long-Term Ecological Research project (PAL-LTER, Ducklow et al. (2012)). For sea ice, we use satellite-based estimates of daily sea-ice concentration with grid resolutions of 25 km (Cavalieri et al., 1996). The datasets used in this study are summarized in Figs. 2-3 and Table 1.

Numerical model
We employ the Massachusetts Institute of Technology general circulation model (MITgcm), which includes dynamic/thermodynamic sea-ice (Losch et al., 2010) and thermodynamic ice shelf (Losch, 2008) capabilities. Following the model configuration from Nakayama et al. (2017), we extract the regional grid from a global LLC270 configuration for the AS and BS regions (Fig. 1).
In the AS and BS domain, horizontal grid spacing is approximately 10 km (Fig. 1). The vertical discretization of the ECCO  Fig. 1) by limiting sea-ice transport between the eastern and central AS. Such a barrier is necessary to simulate 100 sea-ice concentration, and thus air-ocean interaction, closer to observations. The first guess of the model initial state is a simulated 2010 oceanographic condition based on the Green's functions-based solution of Nakayama et al. (2017). Lateral boundary conditions for hydrography, currents, and sea ice are provided by the ECCO LLC270 optimization . The initial guess of surface forcing for 2010-2014 period is from ERA-Interim (Dee et al., 2011). There is no additional freshwater runoff above and beyond the meltwater computed by the MITgcm 105 ice shelf package. The model parameters used for this state estimate are shown in Table 2.
The MITgcm adjoint assimilation system iteratively minimizes a scalar cost function, defined as the weighted least-squares difference between simulation and observations and between prior and adjusted control parameters (e.g., Wunsch et al. (2009);Wunsch and Heimbach (2013);Forget et al. (2015)). The observation weights are spatially homogeneous but depth-varying and defined as the inverse of the simulated variance for potential temperature and salinity at each depth. For example, the 110 estimated error of potential temperature varies from 0.73 • C at the surface to 0.16 • C at a depth of 1000 m. The control vector consists of initial potential temperature and salinity conditions, vertical diffusivity, and time-evolving atmospheric surface boundary conditions (air temperature, specific humidity, precipitation, shortwave radiation, longwave radiation, and eastward and northward winds). Weights for initial temperature and salinity conditions are prescribed to be the inverse variance of the baseline ECCO LLC270 simulation Zhang et al. (2018). The weight for vertical diffusivity is the squared-inverse of the 115 prior value, 5.0×10 −6 m 2 s −1 . Weights for surface boundary conditions are from Chaudhuri et al. (2013). The gradient of the cost function is obtained by integrating the adjoint of the tangent linear model backward in time (Le Dimet and Talagrand, 1986) and is used with the quasi-Newton M1QN3 conjugate-gradient algorithm (Gilbert and Lemaréchal, 1989) to adjust the control variables so as to iteratively reduce the cost function toward its minimum. Despite successful adjoint simulations with particular versions of the sea-ice model (e.g., Fenty and Heimbach (2013a)), the sea-ice adjoint is not used in this study due 120 to persistant instability issues. Sea-ice concentration is instead constrained using a pseudo sea-ice adjoint as is done in Forget et al. (2015). Where the model has an excess (deficiency) of sea ice, extra heat is added to (removed from) the system to bring the sea surface to above (below) the freezing temperature. However, these heat fluxes are only applied when the model has sea ice and observations do not, or vice versa. In this scheme, simulated sea-ice concentration can not be directly optimized. We also note that background horizontal viscosity has to be artificially increased at the early stage of the optimization for model 125 stability and we manually lowered the values of viscosity at iterations 10, 15, and 20 (Tables 2 and 3).
For the static ice shelf component, the freezing/melting process in the sub-ice-shelf cavity is parameterized by the threeequation thermodynamics of Hellmer and Olbers (1989); Jenkins (1991). We use constant turbulent heat and salt exchange coefficients for individual ice shelves, which are already adjusted in Nakayama et al. (2017). However, only for Pine Island and Thwaites, we further modify these coefficients for simulations after iteration 11 (Table 3), as ice shelf melt rates of Pine 130 Island and Thwaites become too large. Changes of these coefficients do not highly alter on-shelf circulation (see Fig. S18 in Nakayama et al. (2018)) and adjustments of these coefficients in addition to optimizations based on adjoint sensitivities is possible.

135
As we initialize the unoptimized simulation (iteration 0) with simulated oceanographic conditions based on Green's functions approach (Nakayama et al., 2017), its 2010 simulated vertical section shows a good agreement with observations ( Fig. 4).  Nakayama et al. (2013)). Similar to Nakayama et al. (2017), slight differences can still be found for WW properties 140 close to the surface (salinity being still too saline (∼0.1 psu)) and PIIS front mCDW properties (∼0.1 psu for salinity and ∼0.2 • C for potential temperature, Fig. 4).
We find, however, that the time evolution of iteration 0 between 2010-2014 does not agree well with observations. For example, oceanographic conditions at the PIIS front in iteration-0 simulation becomes too cold and fresh by ∼2 • C and ∼0.25 psu compared to observations, respectively (Figs. 4 and 5). This is clearly different from observations because WW becomes dense, 145 convects to the bottom, and prevents mCDW intrusions into the PIIS cavity in iteration 0 (Fig. 5). Furthermore, the horizontal section of potential temperature at 552 m depth illustrates the formation of cold and fresh water masses (∼ −1ºC and 34.4 psu) in the vicinity of the PIIS (the red arrows in Fig. 6) in contrast to observations (Fig. 3). This water spreads along the coast and induces unrealistic cooling in the large area of the AS (Figs.6 and 7). Simulated time series of potential temperature and salinity at the PIIS front mooring (Fig. 8) shows that these changes occur as a result of intense cooling by the atmosphere in 150 the austral winter of 2013 and this cooling leads to the reduction of the PIIS melt rate by ∼100 Gt yr −1 (Fig. 9a).

3.2 Model-Observation Differences and improvements
As a result of the iterative optimization, we are able to reduce the cost, which is defined as a sum of weighted model-data difference, by 65% by adjusting initial ocean temperature and salinity, atmospheric surface parameters, and vertical diffusivity ( Fig. 2). The cost reduction occurs quicker in the first 10 iterations (Fig. 2). Throughout the optimization, sea ice and ocean 155 costs are reduced by 59% and 70%, respectively.

Sea ice
In the iteration-20 simulation, spatial patterns of sea-ice concentrations show better agreement with observations. For September, simulated sea ice area over the entire model domain in iteration 0 is larger than observations by 0.08 million km 2 (3.5% difference), while in iteration 20 it is larger by only 0.03 million km 2 (1.2% difference). September simulated sea-ice concen-160 tration is overestimated in iteration 0 at the northern model boundary but becomes much closer to observations in iteration 20 ( Fig. 10). For March in the AS, simulated sea ice in iteration 0 is larger than observations by 0.12 million km 2 (57% difference), and in iteration 20 it is larger by 0.08 million km 2 (37% difference). In the BS, simulated sea ice in iteration 0 is larger than observations by 0.13 million km 2 (144% difference), and in iteration 20 it is larger by 0.05 million km 2 (56% difference).

165
For the AS, there are two major improvements for the oceanographic condition: (1) representation of mCDW intrusions towards the ice shelf cavities and (2) properties of WW. For the BS, we do not include enough observational data in the current version of the ocean state estimate and are not able to judge the capability of our state estimation.  (Jacobs et al., 2011;Nakayama et al., 2013;Dutrieux et al., 2014;Nakayama et al., 2018Nakayama et al., , 2019. The 552-m potential temperature and salinity difference between iteration-0 and iteration-20 simulations are ∼0.5 • C and ∼0.1 psu along the coast of the AS, respectively At shallower depths, the thermocline is located deeper by ∼150 m in the AS in 2014 than in 2010 (Figs.4 and 5), which seems to be consistent with observations (Jacobs et al., 2011;Nakayama et al., 2013;Dutrieux et al., 2014;Heywood et al., 2016). The 222-m salinity difference between iteration-0 and iteration-20 simulations shows freshening by 0.05-0.1 mostly 180 everywhere in the AS (the red arrows in Figs. 7d and e). This is a good improvement as WW tends to become too saline in most numerical models (e.g., St-Laurent et al., 2015;Nakayama et al., 2018) and good representations of surface hydrographic conditions as well as stratification are necessary for the better representation of interannual variabilities. Heat and salt transfer coefficients are kept constant and we do not allow them to change over time for each model iteration. 185 Thus, time series of ice shelf melt rates simply reflect changes of oceanographic conditions in the ice shelf cavities. We note, however, that these coefficients are adjusted once at iteration 11 for Pine Island and Thwaites ice shelves (Tables 3 and 4). In iteration 0, simulated time series of Pine Island and Thwaites ice shelf melt rates show a reduction of ∼100 Gt yr −1 between 2010-2014 ( Fig. 9). In the iteration-20 simulation, however, both time series of Pine Island and Thwaites ice shelf melt rates become rather stable at ∼110 Gt yr −1 but show slight decreasing trends of ∼4 Gt yr −2 and ∼3 Gt yr −2 , respectively. Simulated In general, heat and salt transfer coefficients are already adjusted in Nakayama et al. (2017), and melt rates of ice shelves in the AS and BS are consistent with satellite-based estimates ( Table 4). The interannual variability of the simulated ice shelf melt rates may be too weakened compared to observations possibly because (1) coarse horizontal and vertical grids used in this configuration may not allow the realistic representation of ocean cavity circulation and/or (2) observed reduction of ice 205 shelf melt rates are caused by changes in ice shelf geometry and it can not be simulated in static ice shelf cavity configuration.
Satellite based estimates of time-evolving ice shelf melt rates are required for further comparison.

Sensitivity studies
To investigate the reason for the improvements, we conducted 3 sensitivity experiments (Table 5) as air temperature, precipi-210 tation, and wind are considered to be main drivers of oceanographic variabilities at the PIIS front region. For the NoWindAdj, NoPrepAdj, and NoAtempAdj cases, we re-ran the iteration-20 simulation but excluded adjustments for wind, precipitation, and air temperature, respectively. The total costs are 2.9×10 6 , 3.1×10 6 , 2.9×10 6 , and 4.1×10 6 for iteration 20, NoWindAdj, NoPrepAdj, and NoAtempAdj cases, respectively showing that adjustments of wind and atmospheric temperature play important roles for reducing both sea ice and ocean costs. Sea-ice costs are 1.7×10 6 , 1.6×10 6 , 1.7×10 6 , and 2.9×10 6 for iteration 20, NoWindAdj, NoPrepAdj, and NoAtempAdj cases, respectively . Ocean costs are 1.1×10 6 , 1.5×10 6 , 1.1×10 6 , and 1.2×10 6 for iteration 20, NoWindAdj, NoPrepAdj, and NoAtempAdj cases, respectively. Cost function increases of these sensitivity experiments compared to CTRL are summarized in Table 5, showing that adjustment of wind has the strongest impact on the ocean, while adjustment of air temperature has the strongest impact on sea ice.
Among these sensitivity experiments, the 2014 January mean potential temperature at 552 m depth shows a similar spatial 220 pattern for all cases for open ocean and in the BS (not shown) and differences can only be found in the AS especially at the PIIS front (Fig. 11). Spatially averaged 552 m potential temperatures at the PIIS front (averaged for the region enclosed by red box in Fig 11a) are 0.61 • C, 0.23 • C, 0.56 • C, and 0.53 • C for iteration 20 NoWindAdj" NoPrepAdj, and NoAtempAdj cases, respectively. Vertically integrated heat contents, which are strongly controlled by thermocline depth (Nakayama et al., 2018), reduced by 11%, 5%, and 12%, respectively, for NoWindAdj, NoPrepAdj, and NoAtemAdj cases compared to iteration-20 225 solution. This implies that (1) PIIS front mCDW temperature and thus mCDW pathways as well as strength of intrusions are dominantly controlled by wind and (2) the PIIS front thermocline depth is influenced rather equally by wind, precipitation, and air temperature.

Seasonal and interannual variability
Mooring observations at the PIIS front were conducted from 2009-2014, which provide us potential temperature measurements 230 at various depths (Webber et al., 2017). At depths below 800 m, the observed potential temperature remains rather stable at ∼1ºC (Fig. 12c). At 600-700 m depths, the potential temperature also remains stable at 1 • C and shows gradual cooling and warming between 2010-2012 and 2013-2015, respectively (Fig. 12). Between 2012-2013, however, potential temperature time series shows sporadic emergence of cold watermass (∼-0.5 • C). At 400-500 m depths, the time series of potential temperature fluctuates between -1.8 • C and 0 • C and seasonal and interannual variabilities are large.

235
For iteration 0, we find three major differences with respect to the observations (Fig. 12); (1) simulated potential temperature shows rapid cooling between 2013-2015 and potential temperature at all depths changes from ∼ 1 • C to ∼-1 • C, (2) simulated time series of potential temperature shows sudden emergence of cold water at all depths throughout the simulated period, while it occurs only for shallower depths for observations, and (3) timing of cooling and warming do not agree with observations. This sporadic coolings are likely associated with strong wind events which lead to the formation of cold and dense water and 240 deep convection.
For the iteration-20 simulation, simulated time series show improvements at all depths (Fig. 12). At greater depths (800-900 m), the potential temperature remains rather stable at ∼1ºC consistent observations (Fig. 12c) and the sudden emergence of cold water only occurs at shallower depths. However, the long term and short term variabilities have large differences between observations and the iteration-20 simulation, and the timing of cooling and warming still do not agree with observations. Such 245 differences are also presented in Taylor diagrams based on simulated and observed time series of 400-m and 900-m potential temperature at the PIIS front (Fig. 13). The root mean square differences reduce by 23% and 80% for iterations 0 and 20, respectively. However, standard deviations and correlation coefficients of both the iteration-0 and iteration-20 solutions retain large differences compared to observations. This means that current states of the optimized solution achieve better agreement 8 in terms of mean states but it remains difficult to capture shorter time-scale variability for both time series at the middle and 250 close to the bottom of the water column.
One possible reason for the difference is deficiencies in the simulated sea-ice concentration near the coast. In our current configuration using a pseudo-sea-ice adjoint sensitivity, we are not able to directly adjust sea-ice concentration: simulated sea-ice concentration at the PIIS front remains almost the same between the iteration 0 and iteration-20 simulations (Fig. 12).
More accurate representations of sea ice likely allow us to capture, for example, wind-driven transports of sea ice away from 255 the coast, surface ocean cooling, and sea-ice formation. Such processes likely change the local stratification, which possibly impacts warm mCDW intrusions into the PIIS cavity.

Conclusions
In the previous work, Nakayama et al. (2017) Table 4. Letters E, C, and W denote the submarine glacial troughs located on the eastern AS continental shelf. Transparent white patches (see the red arrow between Do and Cr) indicates the location of grounded icebergs and landfast ice. This white region is treated as a barrier in the sea ice model and we do not allow sea ice exchange crossing this region.The thick black line represents the vertical section shown in Figures 3 and 4. The red circle indicate BSR5/iSTAR9 mooring locations. Both moorings are deployed at the same location off the PIIS.                Table 3. Adjustments to model parameters in addition to optimization using adjoint sensitivities.

Iterations Adjustment
Iteration 10 change background horizontal viscosity from 1000 m 2 s −1 to 500 m 2 s −1 Iteration 11 heat and salt transfer coefficients for PIIS and Thwaites ice shelf reduced by 70% and 63%, respectively Iteration 15 change background horizontal viscosity from 500 m 2 s −1 to 100 m 2 s −1 Iteration 20 change background horizontal viscosity from 100 m 2 s −1 to 20 m 2 s −1 Table 4. Satellite-based estimates of basal melt rate (Rignot et al., 2013) and model mean basal melt rates (2010-2014) for West Antarctic ice shelves for iteration 20. The values of heat transfer coefficient γT used for the optimized simulation are also shown. We use constant turbulent heat and salt exchange coefficients for individual ice shelves, which are already adjusted in Nakayama et al. (2017). However, only for Pine Island and Thwaites (bold), we further modify these coefficients for simulations after iteration 11.