Impact of errors in short wave radiation and its attenuation on modeled upper ocean heat content

Abstract. Photosynthetically available radiation (PAR) and its attenuation with the depth represent a forcing (source) term in the governing equation for the temperature in the oceanic dynamical models. PAR usually comes from the atmospheric model predictions, whereas PAR’s attenuation schemes are internally prescribed (estimated) inside the oceanic dynamical model. We perform sensitivity analyses to investigate the impact that errors in model surface PAR and vertical attenuation of PAR have on the upper ocean model heat content. In the Monterey Bay area, we show that with a decrease in water clarity, the relative error in surface PAR introduces a larger error in the modeled upper 25 m ocean heat content than the same magnitude relative error in the attenuation coefficient. For Jerlov’s type “IA” water (attenuation coefficient is 0.049  m−1), the relative error in surface PAR introduces an error twice as large into the model heat content as the same magnitude relative error in the attenuation coefficient. For the more turbid water Jerlov’s type “III” (attenuation coefficient is 0.127  m−1), the relative error in surface PAR introduces error seven times as large into the model heat content than the same magnitude relative error in the attenuation coefficient. We present how the upper ocean heat content sensitivities to errors in PAR and its attenuation change in space and time. While the sensitivities to the errors in surface PAR are all positive, sensitivities to the errors in attenuation coefficient have positive and negative values, depending on location. They are positive in shallower water for locations on the shelf in the northern part of the bay and negative in deeper offshore waters. Sensitivities derived provide a capability to understand and control the impact of errors in PAR and its attenuation on the upper ocean model heat content predictions.


Introduction
Short wave radiation (SWR) and its attenuation with depth have a major impact on the vertical distribution of oceanic water temperature, dynamical processes, and ocean-atmosphere interactions. SWR is a sum of the ultraviolet-visible wavelengths [350 to 700 nm, photosynthetically available radiation (PAR)] and infrared (IR) wavelengths (700 to 2500 nm). The IR is usually absorbed in the upper few meters of water. However, PAR can penetrate below 50 m in clear water. In existing physical oceanic models, 1-7 the temperature is modeled as a tracer, when its temporal change is described by the advection-diffusion-source equation. The PAR and its attenuation with depth represent a source (forcing) term in the governing equation for the temperature. As a result, in physical ocean models, there are three types of errors related to PAR and its attenuation with the depth: (1) errors in specification of the surface values of PAR, (2) errors in accuracy of the optical model for the PAR attenuation with the depth, and (3) errors in attenuation coefficients of PAR. The remote sensing community has been developing satellite-derived products of PAR and inherent optical properties (IOPs, such as absorption and backscattering, diffuse attenuation coefficient at 490 nm), which are used to estimate PAR attenuation coefficients. [8][9][10] Also, there are well-established optical models, such as HydroLight and EcoLight, [11][12][13] to provide accurate attenuation of PAR with the depth. All of these products can certainly improve modeling of PAR and its attenuation with depth in physical oceanic models. However, at present, there is a very limited use of remotely sensed PAR-related products in oceanic models.
In numerical modeling of oceanic processes, the PAR usually comes from the atmospheric model predictions, and in numerous studies it has been reported that atmospheric models show a tendency to overestimate the SWR due to the underestimation of predicted low-level clouds. [14][15][16][17][18][19] The use of PAR from the atmospheric model is mostly dictated by already in place coupling of the oceanic model to the atmospheric model outputs and also by the fact that the atmospheric model provides forecasts of PAR values that can be used to estimate the corresponding forecast of oceanic conditions, whereas remotely sensed PAR products provide only hindcast or nowcast of PAR values.
The application of optical models based on radiative transfer theory has demonstrated promising results in the idealized channel flow model. 11 However, optical models based on radiative transfer theory have not found wide applications in the numerical oceanic modeling community, mostly because of a high increase in computational cost when they are used with high resolution, data assimilative physical models. For this reason, most existing oceanic physical models 1-7 use a simplified optical model based on an exponential attenuation of PAR with the depth. Specifications of the PAR attenuation coefficient in exponential attenuation schemes are mostly based on: the Jerlov's water-types classification, 20 climatological estimates of attenuation coefficients, or from the biological model predictions of light-absorbing and scattering water constituents. 9,21-25 All of these PAR attenuation coefficient specifications are prone for introducing errors in the attenuation of SWR with the depth.
It is clear that remotely sensed PAR-related products should find wider applications in the oceanic physical models. One way is to assimilate PAR-related products into the oceanic model, when remotely sensed products are merged together with the oceanic model predictions of PAR and its attenuation. In existing data assimilation approaches, this merge is based on estimates of errors in remotely sensed products and corresponding errors in the PAR values and its attenuation in the oceanic model. To achieve this merge, remote sensing and numerical oceanic modeling communities have to better understand the sensitivity of the model heat predictions to the errors in PAR, optical model, and PAR attenuation coefficient.
The sensitivity of the heat modeling to errors introduced by the use of a simplified exponential attenuation scheme was investigated in Ref. 11. In the present paper, our objective is to study the sensitivity of the modeled upper ocean heat content to other mentioned above two types of errors: those related to the specification of the magnitude of PAR at the surface (from atmospheric model) and those related to the incorrect specification of the PAR attenuation coefficient with depth (prescribed in the oceanic model) in the simplified exponential attenuation model.
Among the questions addressed in this study are: What errors (in surface PAR or in its attenuation), and under what conditions, have a greater impact on the modeled oceanic heat content of the area of interest? Under what conditions do the two considered types of errors combine to reinforce the heat content error? Under what conditions do the PAR errors work in opposition, such that they reduce the total error in the heat content? How do the heat content sensitivities to errors in PAR and its attenuation change in time and space?
To address the questions above, we have developed an approach for estimating the sensitivity of the oceanic model heat content to changes in surface values of PAR and its attenuation with the depth. The approach is based on using the adjoint to the oceanic model governing temperature equation. These sensitivities are valid for any model and modeling area where the temperature change is described as tracer with advection-diffusion-source model. Based on derived sensitivities, we have estimated errors in the upper ocean heat content due to errors in surface PAR values and due to errors in the PAR attenuation coefficient.
The structure of the paper is as follows: Sec. 2 describes the modeling system, observations, and the approach for estimating the model heat content sensitivities. Section 3 is devoted to the results. Discussions and conclusions are presented in Sec. 4.

Modeling System
In the present study, the Monterey Bay model [the domain is shown in Fig. 1(a)] is based on the Navy Coastal Ocean Model (NCOM) 2,3 and is triply nested inside the global and regional (California Current) NCOM-based models. 17,18 The model, which is set up on a curvilinear orthogonal grid with resolution ranging from 1 to 4 km, uses a sigma vertical coordinate system with 30 levels. On open boundaries, the model is one-way coupled to a larger scale regional (California Current) NCOM-based model. 17 The Navy Coupled Ocean Data Assimilation system 26 is used for assimilation of observations into the Monterey Bay model. 27 The oceanic model is forced with surface fluxes from the Coupled Ocean-Atmosphere Mesoscale Prediction System (COAMPS) at 3 km horizontal resolution. 28 The COMAPS shortwave calculation considers absorption by water vapor and ozone. Two types of clouds are considered in the radiative transfer calculation: stratiform and cumulus. The Monterey Bay model outputs were validated against temperature, salinity, and currents profiles at moorings, high frequency (HF) radar-derived currents, and aircraft sea surface temperatures (SSTs). 18,27

Approach for Estimation of the Model Heat Content Sensitivities
The governing equation for temperature has the following form in oceanic models: 1,2-7 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 4 8 2 where vector U is the three-dimensional ðx; y; zÞ velocity, k represents the horizontal and vertical diffusivities, Qsr is the SWR at surface, and γ is the function describing the shortwave extinction. γ is 1 at surface and 0 at bottom, ρ 0 is the density, and c p is the specific heat of water. Coordinate z is positive into the ocean and negative up in Eq. (1). Because γ is decreasing in the ocean, the solar radiation in the source term [the last term on the right side of Eq. (1)] is heating the water.
Initial condition for Eq. (1) at t ¼ t 0 has the form E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 6 7 5 where D is the modeling domain and τðx; y; zÞ ∈ D.
In Eq. (1), the Qsr is a sum of the visible-ultraviolet bands (350 to 700 nm, PAR) and IR bands (700 to 2500 nm) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 6 0 9 where f represents a fraction of the visible-ultraviolet bands (PAR) in total solar radiation. The vertical distribution of the PAR and IR is described by two separate exponential functions (see for example Ref. 9) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 5 4 2 where K PAR and K IR are the attenuation coefficients of PAR and IR, respectively. Note that the IR is usually absorbed in the upper meters of water; however, the PAR can penetrate below 50 m in clear water.
Let us consider the heat content J at verification time t f > t 0 in the subdomain Ω (area of interest) of the modeling domain D (our interest is in the upper ocean model heat content) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 4 5 1 where dτ is the element of volume in the integral [Eq. (5)]. As emphasized in Refs. 9 and 29, K PAR varies with depth and horizontally. In Refs. 21 and 22, K PAR is empirically modeled as functions of chlorophyll while in Ref. 9, K PAR is empirically modeled as functions of IOPs (absorption and backscattering at 490 nm), which are derived from the satellite remotely sensed reflectances. Profiles of absorption and backscattering at 490 nm predicted by the biochemical model are used to estimate K PAR in Ref. 25. It is known that there is high level of uncertainty in estimation of chlorophyll and IOPs from remotely sensed reflectances, as well as, based on IOP derivations from the biochemical model simulations. For example, according to Ref. 30, "the satellite data product accuracy generally accepted by the international missions is AE5% for water-leaving radiances and AE35% for chlorophyll a in the open ocean." These relatively large uncertainties make the analytical analysis of errors in the heat content due to errors in K PAR challenging, particularly when K PAR varies horizontally as well as vertically. To simplify the analysis of the heat content sensitivity to errors in K PAR , we assume that K PAR can vary only horizontally. This includes the case when values of f, K PAR , and K IR in Eq. (4) are chosen based on the classification of Jerlov's 20 (as interpreted by Ref. 31), which is currently employed in many oceanic models. 1,3,5,6 Oceanic waters are classified based on its clarity, and five values of K PAR are defined ranging from clear to increasingly turbid water ( Table 1). The case of horizontally varying K PAR also includes a specification of K PAR based on satellite observations of chlorophyll (see for example Ref. 23). Appendix presents the derivations of the errors in the model heat content [Eq. (5)] due to errors in surface PAR and its attenuation with the depth. Equations (15) and (16) provide a representation of errors δJ Q PAR and δJ K PAR in the heat content [Eq. (5)] due to normalized errors ε Q PAR [Eq. (19)] and ε K PAR [Eq. (20)] in surface PAR values (Q PAR ) and in the attenuation coefficient K PAR correspondingly. Equations (15)- (20) provide an estimate of errors in the heat content J [Eq. (5)] at verification time t ¼ t f due to relative errors in Q PAR and K PAR at any time t between initial time t 0 and verification time t f . Also, Eqs.

Satellite ocean color data
The moderate-resolution imaging spectroradiometer (MODIS)-Aqua satellite imagery was processed using the Naval Research Laboratory (NRL) automated processing system (APS). APS is a complete end-to-end system that includes sensor calibration, atmospheric correction (with near-IR correction for coastal waters), and bio-optical inversion. APS incorporates, and is consistent with, the latest NASA MODIS code (SeaDAS). 32,33 The diffuse attenuation coefficient Kd Both above K PAR distributions (August 15 and monthly mean) were interpolated to the model grid.

Mooring data
Observations of winds and PAR from the Monterey Bay Aquarium Research Institute surface moorings M1 (122.02°W, 36.74°N) and M2 (122.40°W, 36.67°N) are used in this study [ Fig. 1(a)]. Near-surface 3-m wind speed and direction were measured by a MetSys monitor. The observed PAR was measured by the spectroradiometer mounted on moorings ∼3 m above the water surface. 35

HF radar surface currents
Surface current observations used in this study were derived from the California Coastal Ocean Current Mapping Program's HF radar network. Surface currents were estimated based on inputs from five HF radar sites for August 2003. Vector currents were estimated on a Cartesian grid with a horizontal resolution of 3 km by computing the best-fit vector velocity components using all radial velocity observations within a radius of 3 km for each grid point each hour. 36

Results
Comparisons of COAMPS predictions with observations from moorings and aircraft surveys were reported in Refs. 18, 27, 28, 37, and 38. The comparisons during August to September of 2003 demonstrated high complex correlations between the observed and COAMPS-predicted wind velocities, and good agreement in the magnitudes and the extent of each observed upwelling event. We estimated the magnitude of the complex correlation coefficient and angular displacement 39 between observed and COAMPS model 10 m wind velocities at the location of the mooring M1 [ Fig. 1(a)]. The estimated value of complex correlation between observed and model wind velocities is 0.79, which is significant at a 95% confidence level. The angular displacement gives the average counterclockwise angle of the COAMPS wind velocity with respect to the observed wind velocity. The angular displacement between model and observed wind velocities is small and equals −11 deg.
At the same time, comparisons of the observed and model-predicted daily-averaged values of PAR [ Fig. 1(c)] show an overestimation of the PAR in the COAMPS predictions, especially during the wind relaxation events of August 20 to 22 and September 1 to 3 (overestimation of mean daily PAR values ranges from 20 to 80 W∕m 2 ). The excessive PAR is likely related to the modeling of low-level clouds. Prediction of extensive low-level clouds in the Monterey Bay area during summer, which is believed to be underestimated by COAMPS during this time period, 28 is a very challenging problem.
To address the research questions raised at the end of Sec. 1, we need to establish "true" values of PAR (Q PAR ) and K PAR against which we estimate errors in the modeled PAR and the model-used values of the PAR attenuation coefficient. The observed PAR values at moorings M1 and M2 can be used to establish magnitudes of errors between the observed and the COAMPS model-predicted values of PAR. In the present study, only satellite-derived values of K PAR are available to us in order to estimate errors in the Monterey Bay model specification of K PAR (from the Jerlov's water-types classification or climatological estimates of attenuation coefficients). This introduces challenges, because satellite-derived K PAR values are prone to their own errors, and satellite images of K PAR are available mostly during clear days, when errors in the COAMPS model predictions of PAR represent the lower bound of possible errors, whereas during cloudy days, when the model predictions of PAR represent the upper bound of errors, the satellite images of K PAR are not available. To mitigate this, we present results for errors in PAR corresponding to the clear day (specifically August 15) and to the monthly mean of PAR for the month of August. In this case, our analysis of the heat content errors will be for the case of the lower bound error in the COAMPS model-predicted PAR values and for the monthly mean error in the modeled PAR. Two model runs were conducted over the considered time period. Run 1 is the model run forced with SWR from the COAMPS predictions. In run 2, the mean daily values of COAMPS SWR were corrected in accord with observed values at the M1 mooring E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 2 6   For both runs, values of f, K PAR , and K IR in Eq. (4) are chosen based on Jerlov's water type "IA" (Table 1). In both runs, no observations were assimilated.
Comparisons of the observed and model-predicted temperatures at 1-and 5-m depths are shown in Figs. 1(d) and 1(e). There is a good agreement between model and observations during strong upwelling favorable winds, when the model dynamics are mostly wind-driven and the solar heating was weaker than other days in August. However, during short wind relaxation events of August 20 to 22 and September 1 to 3, there are substantial differences between model and observed temperatures for run 1, when the model is forced with the surface SWR values without corrections. The run 1 temperatures are much warmer than observed temperatures at both considered depths. However, there is a reasonably good agreement between the model and observed temperatures for run 2, forced with corrected SWR in accord with the observations [Eq. (6)]. According to Table 3, the bias of model temperature predictions at 1-m depth is reduced from 1.83 (run 1) to 0.47°C (run 2), and root-mean-square error (RMS) is reduced at half. This demonstrates that the oceanic model deficiencies in predictions of observed surface and subsurface temperature during wind relaxation events are correlated strongly with the fact that the surface SWR values from the atmospheric model are larger than observed values at mooring locations.
As we have mentioned before, in many oceanic models applications, the K PAR values are chosen either based on the Jerlov classification (specifically, using the attenuation coefficient for the Jerlov type IA in Table 1), or K PAR values are derived from climatological (such as monthly mean) data sets. Table 4 illustrates possible errors in such specifications of K PAR for a specific day-August 15, 2003. According to Table 4, monthly mean values of K PAR represent more turbid, murkier water at both moorings locations in comparison to the K PAR for August 15,2003. At the same time, the Jerlov's type IA water ( Table 4, K PAR ¼ 0.049) represents a more clear, bluish water then for the observed K PAR from the Aqua-MODIS data for the August 15, 2003. Spatial comparisons of monthly mean K PAR with the values for the August 15 are shown in Fig. 2. For most locations in the domain, the monthly mean of K PAR represents more turbid waters than the K PAR for the particular day, August 15, and the error is reaching almost 300% inside the bay.
To    (20)]. Based on those errors and Eqs. (15) and (16), Table 4 shows the ratios of errors δJ K PAR ∕δJ Q PAR in the model heat content. In the first example of   than observed K PAR , i.e., shifted toward more clear water), and the ratio of errors in heat content δJ K PAR ∕δJ Q PAR is also negative (−0.51). In example 2 of Table 5, the error ε Q PAR is positive but larger than in example 1, while the error ε K PAR is negative and is the same as in example 1. As a result, the ratio δJ K PAR ∕δJ Q PAR is still negative (−0.14), and there is also a reduction of the total error in the model heat content.  Table 1. The ratio S K PAR S Q PAR decreases as the water type is changing from clear to increasingly turbid water. This indicates that with the decrease in water clarity, the relative error in surface PAR (ε Q PAR ) introduces a larger error in the model heat content than the relative error in the attenuation coefficient (ε K PAR ) of the same magnitude. For the Jerlov's type 2 water (notation IA in Table 1, K PAR ¼ 0.049), the ratio S K PAR S Q PAR ¼ 0.51, which means that the relative error in the surface PAR introduces error twice as large into the model heat content than the same magnitude relative error in the attenuation coefficient K PAR . For more turbid water Jerlov's type 5 (notation "III" in Table 1 , the relative error in surface PAR introduces error seven times as large into the model heat content than the same magnitude relative error in the attenuation coefficient K PAR . The spatial distribution of the HF radar surface currents [ Fig. 3(a)] and MODIS-Aqua SSTs [ Fig. 3(b)] shows a development of a cold, strong, southward flowing jet along the entrance to the bay during the upwelling event around August 15, 2003. This cold jet separates warmer water masses of cyclonic circulation inside the bay and warmer water masses of the anticyclonic circulation outside the bay [ Fig. 3  and warm anticyclonic offshore eddy, and (3) into the northern part of the bay (upwelling shadow area). This is in accord with circulation during upwelling favorable conditions of August 15, 2003, described above, when the water masses at mooring M1 are mostly influenced by the advection of the upwelled cold water from the upwelling center to the north of the Monterey Bay, interaction between upwelled water of the cold jet and offshore warm anticyclonic circulation, and by the advection of the warmer water from the northern part of the bay along the entrance to the bay. The solution of the adjoint problem formulated in Appendix is nonnegative, and in this case, the sensitivity S Q PAR [Fig. 3(c)] of the heat content [Eq. (5)] to the error in the surface PAR (Q PAR ) is positive in accord with [Eq. (17)]. In contrast to sensitivities to surface PAR errors [ Fig. 3(c)], sensitivities to the errors in K PAR [S K PAR , Figure 3

Discussions and Conclusions
We have presented sensitivity analyses to investigate the impact that errors in model surface PAR and vertical attenuation of PAR have on the upper ocean model heat content. In the Monterey Bay area, we have shown that with a decrease in water clarity, the relative error in surface PAR introduces a larger error in the modeled upper ocean heat content than the same magnitude relative error in the attenuation coefficient. For Jerlov's type IA water (attenuation coefficient is 0.049 m −1 ), the relative error in surface PAR introduces error twice as large into the model heat content as the same magnitude relative error in the attenuation coefficient. For the more turbid water Jerlov's type III (attenuation coefficient is 0.127 m −1 ), the relative error in surface PAR introduces error seven times as large into the model heat content than the same magnitude relative error in the attenuation coefficient. As a result, it is critical to assimilate the available remotely sensed data products of PAR into the atmospheric model (to improve atmospheric model predictions of PAR used by the oceanic model) or to directly assimilate PAR products into the oceanic model by correcting surface heat flux conditions of the oceanic model. In data assimilation, one merges model predictions with observations based on their respective errors. To achieve this merge, remote sensing and numerical oceanic modeling communities have to better understand the sensitivity of the model heat predictions to the errors in PAR and K PAR , which is the paper's objective. The remote sensing community is familiar with errors in PAR and K PAR , however, they are not familiar with how these errors interact in the modeling of the upper ocean heat budget, and how they impact the modeling of the upper ocean heat budget in most existing numerical physical models, and the paper provides this information.
We have presented sensitivities maps of the upper ocean model heat content to the errors in surface PAR and K PAR . While the sensitivities to the errors in surface PAR are all positive, sensitivities to the errors in K PAR have positive and negative values. They are positive in shallower water for locations on the shelf in the northern part of the bay and negative in deeper offshore waters to the west of the M1 mooring. Areas of positive values of both sensitivities indicate that when errors in surface PAR and in K PAR are both positive (or both negative), their impacts on the model heat content error will add up and will result in an increase of the total model heat content error in comparison to the case when errors in PAR and K PAR have opposite signs.
Derived equations for sensitivities are valid for any model where the temperature change is described as a tracer with the advection-diffusion-source model. Based on derived sensitivities, the impact and interaction of errors in PAR and K PAR can be estimated for any modeling domain of interest.
The future of oceanic modeling is in the development of the coupled atmosphere-oceanbiochemical modeling systems with the assimilation of physical, bio-optical observations to improve the initialization of coupled systems and their parameters. The assimilation of available satellite and in situ physical, bio-optical observations should include the assimilation of PAR, as well as K PAR observations. At present, the fidelity of these coupled modeling systems and existing data assimilation approaches are not sufficient for reliable prediction of surface PAR, as well as, of bio-optical properties that impact the modeling of K PAR in the coupled systems (such as estimation of backscattering and partitioned absorption coefficients). With this in mind, it is important to understand how errors in surface PAR (from the atmospheric component of the coupled system), and K PAR (from the biochemical component of the system) will affect the oceanic model predictions, as for example, the upper ocean model heat content. Sensitivities such as those presented here provide a capability to study and understand the impact of different system components on oceanic model predictions.
where S is the surface of the modeling domain D (only water points) and HðsÞ is the water depth at location s ∈ S.
Let us introduce notations δQ PAR ðs; tÞ and δK PAR ðs; tÞ, which are errors in Q PAR and K PAR at time t (t 0 ≤ t ≤ t f ) and at particular location s ∈ S. In this case, the corresponding errors in the heat content J [Eq. (5)] can be written in the following form: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 6 ; 6 7 5 δJ Q PAR ¼ ∂J ∂Q PAR ðs; tÞ · δQ PAR ðs; tÞ; and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 6 ; 6 2 0 δJ K PAR ¼ ∂J ∂K PAR ðs; tÞ · δK PAR ðs; tÞ; where δJ Q PAR and δJ K PAR are errors in the heat content [Eq. (5)] at time t f due to errors in the δQ PAR and δK PAR at time t (t 0 ≤ t ≤ t f ), respectively. Based on Eq. (12), and after simple transformations, we have E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 6 ; 5 4 5 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 6 ; 5 1 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 6 ; 3 5 7 ε K PAR ¼ δK PAR ðs; tÞ K PAR ðs; tÞ ; (20) where ε Q PAR and ε K PAR are the relative errors in Q PAR and K PAR , respectively, while S Q PAR and S K PAR can be interpreted as sensitivities of the heat content J [Eq. (5)] at t ¼ t f to the relative errors in the Q PAR and K PAR . Note that Eqs. (15)- (20) provide estimate of errors (δJ Q PAR and δJ K PAR ) in the heat content J [Eq. (5)] at t ¼ t f due to relative errors in the Q PAR (ε Q PAR ) and K PAR (ε K PAR ) at time t (t 0 ≤ t ≤ t f ), therefore, time between initial time t 0 and verification time t f . In accord with Eqs. (13)- (18), units of errors δJ Q PAR and δJ K PAR , as well as sensitivities S Q PAR and S K PAR are in joules.