Artificial Neural Network Model for Estimating Ocean Heat Content in the Sea Ice-Covered Arctic Regions Using Satellite Data

Ocean heat content (OHC) is an essential parameter to assess Earth’s energy imbalance, global warming, and climate change over the historical record. An accurate estimate of the OHC in the Arctic sea ice regions is challenging due to the lack of in-situ data and satellite-based algorithms. In this study, an artificial neural network-based (ANN) novel approach is presented for estimating OHC changes at various standard depth extents (in line with the World Ocean Atlas-2018 depth levels) in the Arctic sea ice regions based on the relationships of the sea ice thermodynamic parameters from the satellite measurements and in-situ OHC estimates. Because of the potential uncertainty that arises from the inaccessible near-surface oceanic layer in the in-situ OHC estimates, a preliminary ANN model was developed with a set of approximations to account for the in-situ OHC stored in the respective depths within the inaccessible near-surface oceanic layer. The ANN model architecture was optimized for a depth extent of 700 m and adopted for the remaining depths of 20 m, 30 m, 40 m, 50 m, 100 m,150 m, 200 m, 250 m, 300 m, 350 m, 400 m, 450 m, 500 m, 550 m, 600 m, and 650 m. The new model was robust in capturing the spatial, temporal, and depth variabilities of OHC in the sea ice-covered Arctic regions with greater accuracy (mean bias error 0.022 GJ m−2, mean bias percentage error 0.015%, mean absolute error 0.182 GJ m−2, mean absolute percentage error 0.148%, root mean square error 0.24 GJ m−2, R−2 0.94, slope 0.93, and intercept 25.05 GJ m−2). This model is capable of estimating OHC and its temporal trends from satellite data which will have implications for understanding the global climate change and its impacts in the Polar Oceans.


I. INTRODUCTION
The Arctic Ocean is the shallowest of the world oceans, surrounded by many marginal seas, covered by ice throughout the year, and undergoes changes in its sea ice extent with the rise in summer temperatures. Over the recent decades, the Arctic Ocean has received unprecedented global attention as warming temperatures gradually melt sea ice with a declining rate of 10-13% per decade, which makes it the most sensitive region to climate change [1], [2], [3], [4], [5].
The associate editor coordinating the review of this manuscript and approving it for publication was Yongming Li . Persistent sea ice melts in this region have caused changes in the general circulation of the atmosphere and global oceans [6], [7]. This leads to Earth's energy imbalance caused by warming-driven changes such as ice melts and a change in atmospheric water vapor and clouds. The measurements of the ocean heat content (OHC) over time provide the best means of quantifying and analyzing the historical trends in the Earth's energy imbalance [8], [9], [10], [11]. OHC is stored in a slice of the ocean as the absorbed energy and can be used to determine the heat flux in the ocean-atmosphere system and assess climate change. OHC can be estimated as a product of the integrated water column temperature, seawater density, and specific heat capacity from the surface to the desired depth. The amount of heat energy to raise the temperature of a given mass of water varies from absolute zero (0 K) to a particular value (in K); thus, OHC and heat transport calculations prefer the Kelvin scale with a null point of absolute zero for temperature. The OHC calculations use a new variable conservative temperature defined by the TEOS-10 (International Thermodynamic Equation of Seawater-2010) [12], [13] to eliminate the turbulent mixing and pressure effects on the vertical distribution of temperature in the ocean. The OHC calculations use the dynamic density values for each subsurface layer because of its dependence on temperature, salinity, and pressure to quantify ocean stratification with lighter waters near the surface and denser waters in the deep ocean. By considering the above factors, OHC of a given depth layer of the ocean can be defined as the integral of conservative temperature multiplied by the specific heat capacity and seawater density (1), where OHC seawater is given in units of Joules per unit area (J m −2 ), ρ is the seawater density (kg m −3 ), C p is the specific heat capacity (= 3991.867 J kg −1 K −1 ), is the conservative temperature (K) (derived from in-situ temperature, absolute salinity, and pressure), h1 and h2 are upper and lower depths (m), respectively [12], [14]. The Gibbs-SeaWater Oceanographic Toolbox of TEOS-10 is often used to compute the conservative temperature and density for the given in-situ practical salinity, temperature, and pressure values [12], [15].
Traditionally, in-situ-based conductivity, temperature, and depth (CTD) profile measurements are the best means to estimate OHC changes despite some uncertainties due to specific instrument biases, baseline climatology, and spatial data distribution irregularity (mapping method) [16], [17]. However, an accurate reconstruction of regional patterns and long-term global trends in OHC based on in-situ measurements is a challenge owing to the scarcity of data to adequately capture spatio-temporal-depth changes. To minimize the uncertainty associated with sparse observational datasets, ocean reanalysis systems (ORAs) have been developed to provide a complete picture of global ocean variability based on limited observation datasets. The ORAs employ an ocean general circulation model and a data assimilation scheme to synthesize diverse oceanographic products by making use of ocean observations obtained from the in-situ networks and remote sensing systems. Further details of assimilation schemes and approaches to ORA can be found in Stammer et al. [18]. Meyssignac et al. [14] reported a spread of ±0.13 Wm −2 in the OHC trend over 2006-2015 from the surface to 2000 m depth. Higher uncertainty in OHC estimates for the global ocean is plausible for two reasons: i) sparsity of in-situ observational data and undersampled depth profile data that are representative of the global ocean dynamics, and ii) coarseresolution ocean circulation models used in ORAs that do not include the mesoscale eddy processes for the distribution of heat [14], [19]. The ORA estimates of the OHC in the Arctic and its marginal seas diverge quickly with a larger error [20], linked to a dearth or the unavailability of historical in-situ data for assimilation, and thus the poor representation of ocean dynamics [21]. A significant deviation may be caused by limited ice mass balance measurements, near-surface temperature profiles over the ice-covered oceans, and imperfect sampling schemes of ice properties.
Recently, satellite remote sensing has become a promising tool to estimate OHC in space and time, and provides near real-time and wide coverage data that are essentially complementary to observations from in-situ platforms (such as ships, autonomous systems, and fixed stations) for assimilation into numerical models/reanalysis systems to provide a four-dimensional perspective of the global ocean [22], [23], [24]. Numerous satellite-based methods have been developed for accurately estimating OHC changes in the subsurface ocean via empirical statistics, dynamic models, and data assimilation [23], [24], [25], [26], [27], [28], [29], [30]. The potential of these methods is to not only overcome the limitation of data sparsity but to also fill historical data gaps [31]. However, there have been no major studies that have developed and tested remote sensing techniques for estimating OHC in the Arctic region and its marginal seas, due to the lack of in-situ measurement data, the lack of generalization of data, and the difficulty of deriving various ocean parameters and subsurface information in the sea ice-covered Arctic regions [14]. This study is inspired by a critical research gap and presents a novel approach for satellite-based estimates of OHC changes in the sea ice-covered Arctic region. This approach is based on the dynamically linked satellite-derived sea ice thermodynamic parameters with in-situ-based OHC data using an artificial neural network (ANN) model. For this study, in-situ data were obtained through the Woods Hole Oceanographic Institution's Ice Tethered Profiler (WHOI-ITP) program [32], [33] and Cold Regions Research and Engineering Laboratory's Ice Mass Balance (CRREL-IMB) buoy program [34], [35]. The credibility of this approach that links the remote sensing data products with in-situ OHC data using an ANN model, as one of the most popular models of machine learning [31], [36], is fully verified with the test and validation datasets. This ANN model is robust in accurately estimating OHC changes and capturing its spatio-temporal changes with the minimum uncertainty in the determination of decadal trends and regional patterns. The capability of this approach to estimate OHC at 0.25 • spatial resolution from surface to various depth extents is demonstrated using satellite data. This study considers the sea ice-covered ocean north of 60 • latitude that covers the Arctic sea-ice regions, namely the Amundsen basin, Baffin Bay, Barents Sea, Beaufort Sea, Canada basin, Canadian Archipelago, Chukchi Sea, East Siberian Sea, Greenland Sea, Hudson Bay, Kara Sea, Laptev Sea, Makarov basin, and Nansen basin (Fig. 1). In what follows, the database is described in detail (section II). In particular, this includes VOLUME 10, 2022 the calculation of OHC from in-situ profile data and preliminary model simulations, and extraction of local subsets of satellite-derived sea ice thermodynamic data for the insitu matchup locations. Subsequently, the theoretical formulations and sea ice thermodynamics-based ANN model construction and optimization are described for the case of OHC at 700 m depth (section III). The same procedure has been followed for the remaining depth extents. Performance assessment of the OHC model was conducted based on the independent validation dataset (refer to section IV). Further, this work also provides a critical discussion on the potential sources of uncertainties in the OHC estimates (section V).

A. THE IN-SITU OHC DATABASE
The WHOI-ITP is an autonomous instrument designed to provide time-series vertical profile measurements of CTD under the ice at depths between 7 and 750 m. This study utilized the WHOI-ITP Level-3 pressure-bin-averaged CTD profile data with 1 dbar (∼1 m) vertical resolution to calculate the OHC of the profiled water column. For these calculations, profiles with missing or questionable data were excluded to maintain the real trends and variability in the data. The in-situ OHC per unit depth of 1 m subsurface layer was calculated using (1) and integrated over water depth profiled by the WHOI-ITP instrument. The monthly climatological OHC at a spatial resolution of 0.25 • was computed from the objectively analyzed temperature and salinity fields of the World Ocean Atlas-2018 (WOA-18) [37].

1) THE INACCESSIBLE NEAR-SURFACE OCEANIC LAYERS
In this study, the upper ocean column from the surface to 7 m depth was assumed as the inaccessible topmost layer because of the non-availability of the CTD profiling data [20]. Moreover, the upper oceanic column is highly dynamic due to its interactions with the atmosphere. Unlike the open ocean, sea ice (with snow on top of it) essentially forms the boundary between the ocean and the atmosphere for the net heat flux in the ice-covered oceans. This net heat flux at the surface is the external manifestation of changes that occurs from depth-integrated internal heat energy (i.e., OHC). The depth of the ocean column participating in this process may vary in accordance with the intensity of contributing parameters such as solar insolation, surface albedo, outgoing longwave radiation, sea ice freezing and thawing, ocean circulation, and snow cover.
In general, a substance at a temperature above absolute zero on the Kelvin scale must contain heat energy. The heat content accumulated in the inaccessible near-surface oceanic layer comprising snow, sea ice, and a part of the immediate water column underneath must be considered along with the WHOI-ITP sampling layer for establishing a precise relationship between the satellite-based surficial parameters and depth-integrated OHC. To account for the inaccessible nearsurface oceanic layer contribution accurately, a preliminary model was developed to simulate near-surface temperatures and a set of reasonable approximations were made for the other physical oceanographic parameters. This enables an accurate estimate of in-situ OHC in the ice-covered Artic and its marginal seas.

2) A PRELIMINARY MODEL TO SIMULATE THE NEAR-SURFACE TEMPERATURES
The CRREL-IMB buoys provide ice mass balance data for the ice-covered Arctic regions. These buoys are equipped with thermistor strings to measure temperatures from air to sea surface and about 3.8 m depth through the snow or ice and the immediate water column at 10 cm intervals [34], [35]. Using these data, an ANN model has been built that establishes a relationship between in-situ-based nearsurface temperature profile data and satellite-derived surface products such as temperature, sea ice thickness, snow depth, and time of observation (details in section III.B). The Extended Advanced Very High Resolution Radiometer Polar Pathfinder (APP-x) data were used to obtain the satellitebased surface temperature and sea ice thickness products. The Hybrid Coordinate Ocean Model (HYCOM) derived snow depth was obtained from the TOPAZ4 reanalysis system at the Copernicus Marine Service (CMEMS) data portal. Then the ANN model was developed using the data (CRREL-IMB profiles and corresponding satellite data matchups) collected over a time span of 2006-07 and 2010 and validated using the data covering the time period from 2008-09 (refer to Fig. 2). The number of samples used for training (70%) and testing (30%) of this model was 1835. The number of samples used in satellite-based independent validation was 2668 (with a mean bias error of −0.16 K, a mean bias percentage error of −0.05%, a mean absolute error of 1.75 K, a mean absolute percentage error of 0.67%, and a root mean square error of 3.1 K, R 2 0.72).
In addition, the entire dataset was segregated into monthly bins for computing the monthly mean values for the satellitebased input parameters and simulated mean near-surface temperature profiles over the sea ice-covered Arctic region during 2006-2010. In Fig. 3, the simulated near-surface temperature profile data are closely consistent with the CRREL-IMB mean near-surface temperature profile data in the 12 months of a year.

3) ESTIMATION OF IN-SITU OHC IN THE INACCESSIBLE NEAR-SURFACE OCEANIC LAYER
Two possible scenarios were considered with the multiple sub-layers of varying thicknesses (see Fig. 4 & Table 1). The preliminary model was used to simulate temperatures from the surface to 3.8 m depth in both scenarios. Certain approximations and assumptions (italicized in Table 1) were employed for the unattainable temperature measurements in the sub-layers and salinity measurements in the entire inaccessible near-surface oceanic layers. The OHC stored in a given layer of sea ice [38], [39] with a temperature (T i ) was computed using (2), 15 (2)  where T i = in-situ sea ice temperature (K), ρ i = density of sea ice (kg m −3 ), c o = specific heat of fresh ice (2054 J kg −1 K −1 ), L i = latent heat of fusion of ice (3.340 × 10 5 J kg −1 ), and µ = ocean freezing temperature constant (0.054 • C PSU −1 ). The vertical profiles of sea ice salinity (3) were derived from Maykut and Untersteiner [40].
where ω is the normalized coordinate (= depth (z)/ice thickness (h)). This gives the salinity of 0 PSU for the ice-covered surface and 3.2 PSU for the lower boundary layer of the multi-year sea ice. The effect of the sea ice age on salinity variation was neglected when using the above method to derive the sea ice salinity. The sea ice layers with a temperature greater than the melting point (0 • C is considered as melting point) were treated as meltwater or brine depending upon the corresponding salinity (3). The snow layer was approximated as a fresh ice layer for the calculation of OHC snow values.

B. REMOTE SENSING DATA PRODUCTS
Satellite data products used in this study are described in Table 2. The sea ice volume per unit area (SIV) at each grid cell were obtained by multiplying the satellite-derived sea ice concentration with the sea ice thickness data. The sea ice concentration (SIC) is the fraction of a pixel area covered by sea ice and the remaining part of the pixel area was treated as an open ocean fraction. A SIC value of 0 indicates an open ocean area and that of 1 is entirely covered by sea ice. The intermediate values of SIC indicate an open ocean area partially covered by sea ice. The SIC in a particular pixel enables us to understand the level of complexity when analyzing the sea ice thermodynamics within that pixel. In this study, daily SIC data were obtained from the NOAA/National Snow and Ice Data Center (NSIDC) Climate Data Records (CDR). The NSIDC SIC products were estimated using Defence Meteorological Satellite Program's Special Sensor Microwave Imager -Special Sensor Microwave Imager/Sounder (DMSP SSMI-SSMI/S) Passive Microwave data [41], [42].
The APP-x product suite disseminates a wide range of geophysical variables for the Arctic region at two local solar times, 0400 and 1400 hours [43]. Daily sea ice thickness (SIT) and surface temperature (ST) products at 1400 local solar time from the APP-x product suite were used in this study. The sea ice thickness products were estimated using the one-dimensional thermodynamic ice model (OTIM) [44]. The observed mean bias values in OTIM validation during the Scientific Ice Expeditions (1999), Beaufort Gyre Exploration Project (2003)(2004), and at 8 Canadian meteorological stations (2002) are −0.07 m, −0.09 m, and −0.092 m respectively. The OTIM is a physical-statistical hybrid model that uses the parametrizations of residual heat flux and ice physical-thermodynamic processes, whereas the APP-x suite uses a set of linear regression models to retrieve surface temperatures (mean bias of −1.6 K) from AVHRR channels 4 and 5 brightness temperatures and ancillary data such as wind speed and solar zenith angle [43].
Surface and 2 m air temperatures (AT) from satellite observations over the Arctic region were utilized to determine the temperature difference (TD = AT-ST). The air temperatures at 2 m above the ocean surface were obtained from the fifth-generation reanalysis products, generated by the European Centre for Medium-Range Weather Forecasts (ECMWF) system (ERA5). The air temperature products were obtained by interpolating the temperatures between the lowest model level and the surface by considering the atmospheric conditions [45]. Snow depth (SD) data were collected from the TOPAZ4 reanalysis products (mean bias of −0.074 m) at the CMEMS data portal [46], [47], which is based on a recent version of the HYCOM model developed by the University of Miami. The primary dataset containing the in-situ OHC 700 measurements and collocated satellite matchups is divided into two parts: one dataset is used to develop the model and one dataset is to validate the model by ensuring the following requirements.
(i) Spatial distribution: Data points should be welldistributed in the study area to account for all possible geographical variations of OHC resulting from the Atlantic/ Pacific water advection, bathymetry, and river water influx.
(ii) Temporal distribution: Data points should have good temporal distribution so as to render all possible climatic variations caused by solar insolation, wind pattern, precipitation, atmospheric temperature, pressure, and humidity.
The spatial distribution of measurement data points (WHOI-ITP profiles and corresponding satellite data matchups of the sea-ice thermodynamic parameters) covers the entire study region (see Fig. 5). These data points cover a time span of all months from 2006-08 and 2012-15 (used for model development) and 2009-11 (model validation). The total number of samples used for model development (training-70% and testing-30%) is 13932, and for independent validation is 5995. The same criteria is considered for the remaining depth extents as well.

III. DESCRIPTION OF OHC MODEL FOR THE ICE-COVERED ARCTIC OCEAN A. THEORETICAL CONSIDERATIONS
This section presents the theoretical background of the present approach to determine the potential input parameters (satellite-based) to estimate OHC in sea ice-covered Arctic regions by accounting for heat advection by Atlantic and Pacific waters, and heat exchange at the oceanatmosphere (shortwave, longwave, sensible, and latent heat fluxes), ocean-continent (riverine heat flux), and oceanseabed (geothermal heat flux) boundaries in both magnitude and direction (Fig. 6). The sea ice state (thickness, extent, and properties) is sensible enough to OHC changes in the ice-covered Arctic region. A considerable number of studies substantiated that heat at various depths propagates up to VOLUME 10, 2022 FIGURE 5. The spatial distribution of CTD profile data extracted from the WHOI-ITP database. The in-situ measurements used to develop the model are denoted by green circles and those used for model validation by red circles. The ocean regions shallower than 700 m are flagged with a dark grey mask [48].
The heat transfer between the ocean and atmosphere results from the temperature difference between the ocean surface and the air medium, which is greatly controlled by the prevailing SIV at a specific location. A positive value of TD indicates heat gained by the ocean from the atmosphere because of positive (downward) sensible heat fluxes and vice versa. Frozen sea ice leads to heat loss by the ocean (i.e., upward latent heat fluxes from the ocean to the atmosphere), and a similar but opposite condition is observed in the case of thawing. Seasonal variations in radiative and turbulent heat fluxes between the ice-covered ocean surface and the atmosphere are influenced by the amount of sea ice. Vertical heat transports in the ocean are the primary driver of internal heat energy for each grid cell area which has a value of SIC close to unity. In the case of a partial sea ice cover (0 < SIC < 1), significant heat transfer occurs in both vertical and horizontal directions, and the sea ice edges experience large changes in turbulent heat transfer.
The surface area of the sea ice reflects a fraction of incoming solar radiation depending on the amount of albedo (usually within 0.5-0.7 for most sea ice-covered areas) and the level of local insolation (solar irradiance) [61]. The bottom surface of the sea ice provides heat conduction between the ocean and atmosphere through the ice, with significant absorption and transmission of heat energy by the sea ice floes [62], [63]. SIT in a particular grid cell is an essential parameter in calculating the amount of absorption and transmission of heat energy by the sea ice cover. Heat absorption exhibits the phase relation with SIT, whereas an opposite relation exists between heat transmission and SIT. If SIT is equal to or greater than the thermodynamic equilibrium depth (TED) of sea ice, sea ice floes act as an insulating layer by inhibiting heat transfer between the ocean and the atmosphere [61]. It means thermal energy which is available for transfer through a contact surface area is completely absorbed by the sea ice floes.
In addition, the effect of snow cover was considered in the present model to best reproduce the measurement data in the ice-covered Arctic region. Snow has a higher albedo than sea ice (up to a value of 0.9), and it diminishes the quantity of incoming solar radiation reaching the sea ice surface and exhibits similar thermal energy absorption and transmission characteristics as the sea ice layer [61].

B. THE ANN-BASED OHC MODEL
With the theoretical considerations described in the preceding sub-section, an ANN model has been developed to estimate OHC change as a function of the sea ice thermodynamic parameters such as sea ice volume, sea ice fraction, temperature difference, and snow depth. The change in OHC is then combined with monthly climatological OHC at each grid cell area to produce the depth-integrated OHC values from the surface to a stipulated depth (Fig. 7).
The multi-layer perceptron regressor algorithm has been used in this feed-forward ANN model to establish a relation between the input and output parameters [64]. Different combinations of hidden layers count, the number of neurons, feature scaling utilities (StandardScaler, MinMaxScaler, and MaxAbsScaler), activation functions (identity, relu, tanh, and logistic), and weight optimizers (sgd, adam, and lbfgs) were experimented on a depth extent of 700 m (OHC 700 ) to arrive at an optimum ANN architecture with better validation metrics and the same has been adopted for the remaining depth extents. The number of hidden layers used in this architecture is 5 (it is 4 for the preliminary ANN model, refer to section II.A.2), with 100 neurons in each of the hidden layers.
The StandardScaler class of Scikit Learn was employed for feature scaling, and it standardized the data to have a mean of 0 and a variance of 1. For this, it computes the mean and the standard deviation of the training data, and applies the following transformation to each of the samples (both training and testing data).

Standardized or scaled data =
Raw data − Mean Standard deviation (4) In the first iteration (forward step), the input layer holds the scaled data, and this supervised learning model assigns the  random weights (W ij ) for all the connections in the network. The neurons at each ANN layer (except the input layer) collect the previous layers' data and compute the weighted sum using the summation function. Subsequently, it transforms the data using the rectified linear unit (relu) activation function (except the output layer). In the case of the preliminary ANN model, the identity function was outperformed the remaining activation functions.
relu activation function f (x) = maximum(0, x) (5) identity function f (x) = x (6) The computed output was compared against the actual output data at the output layer, and the corresponding squared error was calculated.
Bias or Error = Computed data − Actual data The backpropagation algorithm passes this error from the output layer to the previous layers and employs the Adam optimizer (a stochastic gradient-based optimizer) for weight optimization. The process is repeated until the improvement in error is greater than the specified tolerance. A tolerance value of 0.000001 was considered in this work.
Weight modified = Weight previous − learning rate × ∇Error (8) where the learning rate is the step-size in updating weights = 0.001, and the gradient of the error is computed with respect to the weights in the previous iteration. This ANN model (Fig. 8) took 39 iterations to reach convergence with the root mean squared loss of 0.19 GJ m −2 . Separate ANN models were built for each standard depth level considered, and unified as a generalized model with the help of a python script. This python script loads the concerned individual model based on the bathymetry level at that particular location and estimates the OHC for the given input data of sea ice thermodynamic parameters.

IV. PERFORMANCE ASSESSMENT
The performance of the model is assessed by comparing the model outputs with the measured values of the dependent variable (OHC) using certain global statistical metrics. These statistical metrics include mean bias error (MBE), mean bias percentage error (MBPE), mean absolute error (MAE), mean absolute percentage error (MAPE), root mean square error (RMSE), squared (Pearson) correlation coefficient (R 2 ), slope, and intercept. Some of these metrics are defined below The performance of the model was assessed by comparing the model estimates with the independent in-situ OHC data over a wide range of conditions ( Fig. 9 & Table 3). The observed statistical metrics and scatter plots demonstrating the applicability of sea ice thermodynamics-based ANN model in estimating OHC changes at various depth extents. Further, the model-derived OHC values are compared with the ORA-derived OHC values from the temperature and salinity profile data of the Multi Observation Global Ocean ARMOR3D L4 analysis system (for the standard depths considered) [65]. The 3-year mean data of the period 2009-11 is considered in this comparison (N = 81,802) and a relative difference of −0.24 GJ m −2 (−0.09%) is observed between the model-derived and ORA-derived OHC estimates computed using the python script (refer to section III.B). The observed difference is because of the discrepancies between the various data sources and methods adopted respectively. However, detailed insights on the potential uncertainties and possible solutions to avoid/minimize the error contributions are discussed in the subsequent section.

V. DISCUSSION
The spatio-temporal changes and regional patterns in OHC can be accurately estimated using subsurface temperature measurement data which are collected over decades and cover VOLUME 10, 2022 a widespread area across latitudes and various depth ranges of OHC accumulation. This implies that inadequate measurement data or inaccurate models could introduce large biases in the long-term OHC trends and regional OHC patterns [66]. The in-situ-based OHC estimates are generally associated with an accuracy of ±0.11 Wm −2 , which is within the accepted limit of ±0.3 Wm −2 [14]. A higher uncertainty is generally associated with traditional satellite-based OHC estimates due to the temporal correlation errors and significant constraints of the OHC models due to the sparsity of in-situ data. The remote-sensing techniques often produce inaccurate OHC estimates because of their limitations to account for all of the influencing factors that affect OHC changes and minimize the sources of uncertainty in the icecovered Arctic ocean.
There has always been difficult to obtain subsurface CTD profile data in the ice-covered Arctic ocean. A rough decadal mean estimate of OHC change based on the WOA-13 data showed ∼4% of global OHC change in the Arctic region [67]. The ocean reanalysis system has the limited capability to capture the large-scale variability of OHC change in this region. The limited in-situ observing systems and different strategies to calculate the OHC trends could introduce the different sources of uncertainty (which change on both spatial and temporal scales) over recent periods.
This study elaborates the sources of uncertainty in OHC changes based on our sea-ice thermodynamics model and provides possible solutions to mitigate a statistically significant contribution of error in OHC estimates in the ice-covered oceanic regions. In-situ observations showed a significant variation of OHC on hour timescales caused by the variation of solar insolation, sea surface winds, humidity, cloud cover, precipitation, and other factors. To elucidate this variation of OHC, observations from both in-situ and satellite instruments were assimilated in this study, which used the WHOI Level-3 ITP CTD data at a time interval of about 0000, 0600, 1200, and 1800 hours and the satellite-derived products at 1400 hour. Our analysis showed that a significant part of the uncertainty in the OHC estimates comes from the temporal drift errors under natural ambient conditions. There are two possible solutions to reduce this temporal drift error: i) insitu measurements are to be made within 1 hour of satellite overpass (matchups), and ii) in-situ measurements are to be collected simultaneously for all of the model parameters. Simultaneous CTD and sea-ice mass balance measurements are needed to minimize the temporal drift error in the OHC estimates.
The NOAA/NCEI climatological maps of WOA-18 were produced by implementing the objective analysis interpolation algorithm on the World Ocean Database-2018 (WOD). The Arctic ocean is the least covered basin and has sparse in-situ data compared to the remaining Ocean basins [20]. This sparsity may cause unavoidable biases while executing objective analysis interpolation for the Arctic region. These biases are expected to decrease by updating the WOD with newly collected in-situ data. The NOAA/NCEI is periodically updating the WOD and releasing the updated WOA for the scientific community. Thus the most recent WOA version may be used to improve depth-integrated OHC computations.
The present study employed a set of reasonable approximations and assumptions on the calculation of in-situ OHC in the inaccessible top ocean layer. This led to a small uncertainty in in-situ-OHC estimates in the subsurface layers. However, to reduce this uncertainty, an extensive effort is required to create an in-situ CTD profile database with a large number of samples for different layers such as snow, sea ice, and immediate water column. In addition, prototype modelling studies are needed to better understand sea ice thermodynamics, temperature, and salinity variability in the inaccessible top ocean layer and refine the approximations and assumptions to achieve a higher accuracy level.
The accuracy of the satellite-derived products is critical for developing a reliable model and assessing its accuracy for real-world applications. It is anticipated that future research will improve the accuracy of the sea-ice thermodynamic parameters such as SIC, SIT, SD, ST, and AT. More precisely, there is a lack of satellite-derived snow concentration products to calculate snow volume for this OHC model; a major improvement in the OHC estimates would be achieved by using the snow volume along with SD. These improved satellite-derived products will subsequently enhance the accuracy and efficiency of the present model to estimate OHC changes in the Arctic region.

VI. CONCLUSION
Accurate determination and reconstruction of OHC estimates can be achieved by means of the well-distributed CTD profiles in time (decades and longer) and space (sufficiently widespread) extents. They enable us to detect and analyze the spatio-temporal changes in Earth's energy imbalance, and its consequences on climate change. However, there have been changes evolving with regard to instrumentation, geographic range, and depth coverage which can introduce uncertainty in the determination of long-term OHC trends and its regional patterns [66]. There is also the data sparsity and the lack of common approaches for a reliable OHC estimate over the ice-covered Arctic region, leading to poorly constrained models to capture spatio-temporal variability of OHC in the region. In this study, a sea-ice thermodynamics-based model was developed based on ANNs. This model is robust to the sparse observation data and provides an accurate estimate of OHC change from the surface to various standard depths at a spatial scale of 0.25 • . It accounts for the major sources of uncertainty and minimizes the geographic noise arising from sparse observation data. This model was trained to estimate the depth-integrated OHC using in-situ observation data. The model architecture was optimized to estimate OHC changes from a set of satellite-derived sea-ice thermodynamic data, including sea ice thickness, sea ice concentration, snow depth, surface temperature, ambient air temperatures, and climatological OHC. For this model, CTD profile data were obtained through the WHOI-ITP program to generate an in-situ database for the upper layer from the near-surface (7 m depth) to various standard depth extents. The OHC of the inaccessible top ocean layer was then estimated by a preliminary ANN model, and a set of reasonable approximations and assumptions. This preliminary model simulated the near-surface temperature profiles (from surface to 3.8 m) of the top ocean layer covered by snow, sea ice, and a portion of the immediate water column. The sea ice mass balance measurements and nearsurface temperature profiles obtained from the CRREL-IMB program were used to build this preliminary ANN model. The new OHC model is a promising tool to estimate consistent spatial and temporal OHC changes in the ice-covered Arctic region.
At present, a full-depth estimate of OHC is less feasible due to the sparsity of in-situ data at greater depths and the OHC estimates are associated with a considerable uncertainty for the layers from the surface to 700 m. There are directions to improve the uncertainty in satellite OHC estimates: i) to improve the uncertainty in the in-situ data based on the calibration and cross-validation with multiple platform observations, ii) to optimize the ANN model with more independent validation data to produce a better confidence in the modelderived OHC, iii) to improve the satellite-based estimates of thermodynamic parameters, and iv) to conduct further studies on the heat fluxes which are critical to the thermodynamics based model in the ice-covered top ocean layer.

DATA AVAILABILITY
The Woods Hole Oceanographic Institution's Ice Tethered Profiler program data (Level-3 averaged at 1-dbar vertical resolution) can be downloaded from ftp://ftp.whoi.edu/ whoinet/itpdata. Ice mass-balance measurements and near-surface temperature data are available at the Cold Regions Research and Engineering Laboratory's Ice Mass Balance buoys program website ( [35] and http://imbcrrel-dartmouth.org/archived-data/). Monthly climatological CTD profile data at 0.25 • spatial resolution can be extracted from the World Ocean Atlas-2018 ( [37] and https://accession.nodc.noaa.gov/NCEI-WOA18). Satellitederived (DMSP SSMI-SSMI/S Passive Microwave data) Daily sea ice concentration data at 0.25 • spatial resolution can be downloaded from the NSIDC data portal ( [42] and https://doi.org/10.7265/N59P2ZTG). Daily sea ice thickness SIT and surface temperature ST products at 1400 local solar time can be obtained from the NOAA Climate Data Record of the Extended AVHRR Polar Pathfinder product suite (Key et al. [43] and http://doi.org/10.25921/AE96-0E57).