Leveraging Uncertainty Quantification to Design Ocean Climate Observing Systems

Ocean observations are expensive and difficult to collect. Designing effective ocean observing systems therefore warrants deliberate, quantitative strategies. We leverage adjoint modeling and Hessian uncertainty quantification (UQ) within the ECCO (Estimating the Circulation and Climate of the Ocean) framework to explore a new design strategy for ocean climate observing systems. Within this context, an observing system is optimal if it minimizes uncertainty in a set of investigator‐defined quantities of interest (QoIs), such as oceanic transports or other key climate indices. We show that Hessian UQ unifies three design concepts. (1) An observing system reduces uncertainty in a target QoI most effectively when it is sensitive to the same dynamical controls as the QoI. The dynamical controls are exposed by the Hessian eigenvector patterns of the model‐data misfit function. (2) Orthogonality of the Hessian eigenvectors rigorously accounts for redundancy between distinct members of the observing system. (3) The Hessian eigenvalues determine the overall effectiveness of the observing system, and are controlled by the sensitivity‐to‐noise ratio of the observational assets (analogous to the statistical signal‐to‐noise ratio). We illustrate Hessian UQ and its three underlying concepts in a North Atlantic case study. Sea surface temperature observations inform mainly local air‐sea fluxes. In contrast, subsurface temperature observations reduce uncertainty over basin‐wide scales, and can therefore inform transport QoIs at great distances. This research provides insight into the design of effective observing systems that maximally inform the target QoIs, while being complementary to the existing observational database.

co-design within the oceanographic community. First, it gives insights into the physical mechanisms that govern optimal design strategies; and second, it quantitatively assesses redundancy and optimality of an (existing or future) observing system.
To place our technique into context, we briefly recall existing formal approaches to observing system design.
Observing System [Simulation] Experiments (OS[S]Es; Fujii et al., 2019) are the most common computational tools in oceanography to support observing system design (e.g., Balmaseda et al., 2007;Gasparin et al., 2019;Griffa et al., 2006;Halliwell et al., 2017). OSEs are limited to evaluate existing observing systems, whereas OSSEs can test the skill of proposed future observing systems. The design strategy to be tested in an OSSE has to be specified by the investigator. Once a region is targeted for monitoring, the proposed observing system design (to be tested in the OSSE) is typically guided by best available knowledge of both local hydrographic properties and local dynamical balances (Hirschi et al., 2003;Li et al., 2017;Perez et al., 2011). An example are the Atlantic trans-basin mooring arrays OSNAP (Li et al., 2017), RAPID (Hirschi et al., 2003), and SAMBA (Perez et al., 2011), which target monitoring of the meridional overturning circulation. Key components of each design are western and eastern boundary moorings, which allow geostrophic transport estimates across each trans-basin section. Although these local considerations support practical local monitoring, it is possible that similar constraints could be obtained elsewhere, perhaps with an instrument configuration more easily sustained, at reduced cost, or less susceptible to noise. This opportunity arises from the appreciation that observed variability at any given location is rarely a purely instantaneous response to local forcing. Instead, it is the superposition of phenomena originating in distant regions and at distinct times, communicated through the ocean by advection, diffusion and wave propagation (Heimbach et al., 2011).
Exploring the possibility of remote constraint is essential for truly optimal observing system design. Adjoint models have proven valuable for fully mapping the local and remote origins and pathways of variability in targeted quantities, for example, meridional overturning at given latitudes (Heimbach et al., 2011;Köhl, 2005;Pillar et al., 2016;Smith & Heimbach, 2019). Exploiting the rich information exposed by an adjoint model, a number of adjoint modeling techniques have previously been used to inform ocean observing system design, for example adjoint sensitivity (Heimbach et al., 2011;Masuda et al., 2010), observation sensitivity (Köhl & Stammer, 2004;Moore et al., 2011), and singular vectors (Fujii et al., 2008;Zanna et al., 2012). Despite giving valuable insight into where observations may be useful, none of these latter techniques provide a measure of redundancy versus complementarity, nor of optimality of an observing system. These obstacles are overcome by Hessian UQ: an adjoint-based technique embedded in a variational data assimilation system. The Hessian matrix (composed of second derivatives) of the cost function J captures the curvature of J with respect to the control variables, and allows one to calculate how much uncertainty is reduced with any changes applied to the observing system (Thacker, 1989). In contrast to the previous adjoint modeling techniques named above, Hessian UQ accounts for data redundancy. It also provides a measure of optimality: the more uncertainty an observing system reduces in a defined target quantity (on a scale of 0%-100%), the closer it is considered to being optimal for the defined target.
Hessian UQ has been routinely applied in numerical weather prediction (NWP; Leutbecher, 2003) and, more broadly, in computational science and engineering (CSE; Bui-Thanh et al., 2012), but it has only seen limited use in the oceanographic community. Previous studies have applied Hessian UQ after severely reducing the dimension of the space of uncertain parameters in an ad-hoc manner (Kaminski et al., 2015;Kaminski et al., 2018), or in the dual form of "representers" (Bennett, 1985;Moore et al., 2017;Zhang et al., 2010). These examples have focused on regional ocean settings and on daily to monthly time scales. In this study, we take a step toward fully exploiting Hessian UQ to design ocean observing systems that are targeted at climate monitoring in a global context. To this aim, we apply Hessian UQ within the global ocean state estimation framework of the ECCO consortium , and elucidate oceanic teleconnections that communicate observational constraints over basin-scale distances and monthly to interannual time scales.
In ocean climate research, the goal of an observing system is usually to accurately estimate certain quantities of interest (QoIs): forecasts or climate indices that are difficult or impossible to observe directly. Examples of QoIs include transports across certain oceanographic passages, ocean heat content near the polar ice sheets, regional sea level anomalies, or future sea-ice extent. We therefore focus on the information that an LOOSE AND HEIMBACH 10.1029/2020MS002386 2 of 25 observing system contains about a given QoI, here referred to as the observing system's "proxy potential" for the QoI on a scale of 0%-100% (Loose et al., 2020). Proxy potential is defined by way of Hessian UQ, as the reduction in QoI uncertainty that would be achieved if the observing system was added to the ocean state estimate. Importantly, proxy potential can be assessed not only for existing but also for future observing systems, because it does not require the actual measurement values of the observations (only their locations, times, types, and uncertainties). Loose et al. (2020) provided interpretations of Hessian UQ and proxy potential for idealized cases, in which an observing "system" consists of only a single and noise-free observation. Then, the observation's proxy potential for a QoI reflects the degree to which adjustment mechanisms are shared between the observation and QoI. In this simple case, proxy potential can be understood as the dynamical analog of statistical correlation (squared) between observation and QoI, with the important distinction that proxy potential accounts only for covariability that has dynamical underpinnings. The goal of this study is to leverage Hessian UQ to generalize the notion of proxy potential introduced by Loose et al. (2020) in three important ways (Section 2): first, by extending this concept from a single observational asset to full observing systems; second, by quantifying observational redundancy versus complementarity; and third, by accounting for observational noise.
We illustrate the concepts of Hessian UQ and proxy potential in a North Atlantic case study (Section 3). To provide a clear understanding of Hessian UQ, our case study focuses on observing systems that are comprised of only a few observations. We then discuss how our approach and the dynamical insights obtained generalize to the design of full-fledged observing systems, including thousands to millions of observations (Section 4).

Ocean State Estimation
Ocean state estimation optimally fits an ocean general circulation model (GCM) to the available observations in a dynamically and kinematically consistent way. For this, one solves an inverse problem: given an observing system (gray box, Figure 1), one adjusts the control vector   1 ( , , ) T N u u u (green box, Figure 1), such as to minimize the scalar cost function The control variables u 1 , …, u N (i.e., the elements of the control vector u) are the uncertain input variables of the model, and consist not only of initial conditions (as common in NWP), but also of atmospheric forcing variables and uncertain model parameters (green box, Figure 1). The function J misfit (u) measures the misfit between the vector of actual observations,    Figure 1), given the input variables u. The function J prior (u) penalizes deviations from a first-guess u 0 of uncertain inputs. The M × M matrix R and N × N matrix B are chosen error covariances, spelling out the assumption that observational noise and prior uncertainties follow the Gaussian distributions ( , ) 0 R  and 0 ( , ) u B  , respectively (Tarantola, 2005).
The solution of the inverse problem is the minimizer of the cost function, u min = min u J; that is, a choice of control variables. The ocean state estimate itself is obtained by running the GCM with inputs u min .

Inverse Uncertainty Propagation
To quantify uncertainties in the solution u min of the inverse problem, one propagates observational information and uncertainty along path (UQ1; Figure 1). This inverse uncertainty propagation results in the LOOSE AND HEIMBACH 10.1029/2020MS002386 3 of 25 posterior probability distribution of the control variables, given the observations. In practice, it is not feasible to compute the full posterior probability distribution, nor to map this distribution onto the full ocean state space. We therefore need to appeal to approximation methods.
The posterior probability distribution of the control variables can be approximated by the Gaussian Here, is the set of orthonormal eigenvectors v i with associated nonzero eigenvalues λ 1 ≥…≥ λ M′ > 0 of the misfit Hessian: Starting from an observing system (gray box), inverse uncertainty propagation along path (UQ1) reduces the uncertainty in the control variables (green box), see Section 2.2. A subsequent forward uncertainty propagation along path (UQ2) reduces the uncertainty in a chosen quantity of interest (QoI, purple box), see Section 2.3. Green and black arrows indicate propagation of prior and posterior uncertainty, respectively. The degree to which the observing system reduces uncertainty in the QoI, via a composite uncertainty propagation along paths (UQ1) and (UQ2), is referred to as the observing system's proxy potential for the QoI (Section 2.4).
In Equation 3, the entries of the M × N matrix A are the sensitivities u, and can be thought of as a nondimensionalization if B is diagonal. In summary, Equation 2 phrases the posterior uncertainty P as the prior uncertainty B, reduced by any information {v i , λ i } obtained from the observations. Equation 2 has been known and used in the NWP and CSE communities for many years (see, e.g., Bui-Thanh et al., 2012;Leutbecher, 2003). A self-contained derivation is relegated to the supporting information (Text S1).
Next, we inspect the set {v i , λ i } in more detail as it fully characterizes the information obtained from the observations. The eigenvectors Put differently, the only data-informed direction is spanned by the prior-weighted sensitivity vector  /2 Obs T u B (Equation 5), where "prior-weighting" is through multiplication by B T/2 . Similarly, for an observing system with more than one observation (M > 1), the data-informed subspace of the control space is spanned by the M prior-weighted sensitivity vectors B T/2 ∇ u Obs 1 , … , B T/2 ∇ u Obs M . To obtain the eigenvectors of the misfit Hessian, one has to orthonormalize and rotate these M vectors within the data-informed subspace (Appendix A). In particular, the eigenvectors of our misfit Hessian-which contains the second derivatives of J misfit -are fully determined by first (rather than second) derivatives of the observed quantities, that is, by ∇ u Obs i .
For M = 1, the observational noise, ɛ 2 , appears in the denominator of λ 1 (Equation 5). In particular, for vanishing ɛ 2 , the eigenvalue λ 1 tends to infinity. This fact generalizes to the case M > 1: in the limit of vanishing observational noise (  0 R ), the eigenvalues λ i of the misfit Hessian (Equation 3) tend to infinity, That is,

Forward Uncertainty Propagation
To assess the observational constraints on a QoI, the inverse uncertainty propagation along path (UQ1; Figure 1) has to be followed by a forward uncertainty propagation along path (UQ2; Figure 1). In other words, we quantify how the uncertainty reduction in the controls, due to the new observational information, reduces uncertainty in the QoI, a diagnostic of the model evaluated at u min . Forward propagation of prior uncertainties (B, dotted green arrow) and posterior uncertainties (P, dotted black arrow) along path (UQ2) results in the prior and posterior QoI variances (see Isaac et al., 2015, or Text S2 in the supporting information):  10), orthogonally decomposed into q = q obs + q null . The datainformed component, q obs , is the projection of q onto the data-informed subspace. The component q null lies in the nullspace, that is, the subspace that is not informed by the data. Parts of the unit sphere of the control space are displayed in black, and q has unit length. The larger the radius of the orange dashed circle, defined by the length of q obs , the higher the dynamical proxy potential of the considered observing system for the QoI.
We infer the prior-to-posterior reduction in QoI uncertainty, relative to the prior uncertainty: The second equality in Equation 9 holds by virtue of Equations 8 and 2. Here, are the eigenvectors and eigenvalues of the misfit Hessian (Equation 3), • denotes the "dot" (or Euclidean inner) product between two vectors in N  , and The unit vector q is of key interest: it is the direction within the control space to be constrained in order to inform the QoI. It can be written as the orthogonal decomposition q = q obs + q null (Figure 2c). q obs is the component that lies in the data-informed subspace, given by the projection The component q null lies in the orthogonal complement of the data-informed subspace: the null space, that is, the subspace that is not informed by the data. Uncertainty is only reduced along the data-informed component, q obs , not along the nullspace component, q null .

Dynamical and Effective Proxy Potential
Relative reduction in QoI uncertainty,  2 QoI Δ (Equation 9), rigorously quantifies the dynamical constraints of an observing system (gray box, Figure 1) on a QoI (purple box, Figure 1), as the result of the composite uncertainty propagation along paths (UQ1) and (UQ2). We refer to  2 QoI Δ as the proxy potential of the observing system for the QoI (Loose et al., 2020). Building on Equation 9, we distinguish between dynamical proxy potential and effective proxy potential of the examined observing system for the QoI. Recall that M′ ≤ M is the number of independent data constraints, characterized by the eigenvectors and eigenvalues ). Geometrically, DPP is equal to the squared length of q obs , the data-informed component of q in control space ( Figure 2c). Note that EPP is smaller than DPP, because all factors η i = λ i /(λ i + 1) are smaller than 1. For vanishing observational noise, EPP approaches DPP, since all eigenvalues λ i tend to infinity (Equation 7), and consequently all factors η i = λ i /(λ i + 1) tend to 1 (see also Appendix B). The bounds for DPP and EPP correspond to the cases for which the observing system provides no constraint (EPP = DPP = 0), and for which it serves as a perfect proxy for the QoI, in the case of noise-free observations (DPP = 1) and noisy observations (  EPP 1).

If the observing system under consideration consists of only a single observation (M = 1), Equation 11
simplifies to 

Application to the North Atlantic
We illustrate the concepts of Hessian UQ and proxy potential in a North Atlantic case study. Section 3.1 describes our experimental setup, including our choice of QoI and observations. We then assess proxy potential of the observations for the QoI, for the cases of noise-free observations (DPP, Section 3.2) and noisy observations (EPP, Section 3.3).

Experimental Setup
Our experimental setup coincides with the one described in section 3.1 of Loose et al. (2020) and is embedded in the ECCO version 4, release 2 (ECCOv4r2, Forget et al., 2015) state estimation framework. We use the Massachusetts Institute of Technology GCM (Adcroft et al., 2018;Marshall et al., 1997), in a global configuration, at a nominal horizontal resolution of 1°, and with 50 vertical levels. The linear sensitivities of the QoI and observed quantities to all control variables (Equations 5 and 10, Appendix A) are computed using the respective adjoint models generated through algorithmic differentiation with the commercial tool Transformation of Algorithms in Fortran (TAF; Giering & Kaminski, 2003).
Our QoI is heat transport across the Iceland-Scotland ridge (black line, Figure 3), denoted by HT ISR . We study four different hypothetical temperature observations in the North Atlantic, located inside the yellow dots in Figure 3, and labeled by θ A , θ B , θ C , θ D . Observations θ A and θ C are in the Irminger Sea at (40 °W, 60 °N), observation θ B off the Portuguese coast at (12 °W, 41 °N), and θ D in Denmark Strait at (28 °W, 66 °N). θ A , θ B , and θ D are subsurface observations, situated at 300 m depth, whereas θ C is a surface observation.
Exact definitions for the model calculations of HT ISR and   can be found in Loose et al. (2020). We quantify dynamical and effective proxy potential of the five-year mean of the observations for the five-year mean of our QoI, for zero lag. Sensitivities of the QoI (Equation 10) and observations (Equation 5) to the control variables are computed from the final five years (2007-2011) of the ECCOv4r2 state estimate.
As in Loose et al. (2020), the control variables are comprised of the time-mean, spatially-varying forcing fields Q net,↑ , EPR, τ x , and τ y (Table 1)  . The length of the control vector, N, is 6 (10 )  , equal to 4 times the number of model surface grid cells covering the global ocean (next to last column, Table 1). We assign each of the four forcing fields, F m (i, j), a spatially constant prior standard deviation, ΔF m (last column, Table 1; see also Loose et al., 2020). Moreover, prior cross-correlations are set to zero, corresponding to the assumption that the decorrelation length in the surface forcing is less than the grid scale (∼1°). These choices result in a diagonal N × N prior covariance matrix, B, whose diagonal is populated by the 2 (Δ ) m F assigned in Table 1. The rationale for our simplified choice of control variables and prior covariance matrix is that it allows a more concise presentation of the results, facilitating a clear understanding of the methodology presented in Section 2 and its underlying dynamical and mathematical concepts.

Noise-Free Observations
We begin by exploring the DPP of the temperature observations   , ⋆ = A, B, C, D, for our QoI, HT ISR . DPP is equal to the relative uncertainty reduction in HT ISR that is achieved when adding   to the underlying state estimation framework, in the limit of vanishing observational noise (Equation 11). We first study the DPP of each observation individually (Section 3.2.1), and then the DPP of observing systems that are formed by multiple   (Section 3.2.2).

Degree of Shared Adjustment Mechanisms
The DPP of each observation   for HT ISR quantifies the degree to which the direction of interest, q, projects onto the   -informed subspace of the control space (Equation 11 and Figure 2c). We denote the (one-dimensional)   -informed direction by v  . As an example, Figures 4a-4e show q (Equation 10) and v  (Equation 5), restricted to their τ y (meridional wind stress) components: Here, the normalization factors   HT 6TW B and  B  , whose values are reported in Table 2, are prior uncertainties, computed according to Equations 8 and 6, respectively. In Equation 13, prior-weighting, that is, multiplication by B T/2 , has simplified to an element-wise multiplication of the sensitivity vectors by the constant Δτ y ( Table 1) The Q net , EPR, and τ x vector components are not shown, but their relative contributions to the magnitude of the vectors q and v  are illustrated in Figures 4f-4j. The relative contribution of each forcing field F m to q is computed via and similarly for v  . relative importance of each forcing field, F m , for impacting HT ISR and   , respectively. Wind forcing (τ x , τ y ) is more influential than buoyancy forcing (Q net , EPR) for HT ISR and for all observations   , except for the surface observation θ C , which is highly sensitive to local air-sea heat fluxes (Figure 4i).
The yellow labels in Figure 4 present the DPP (Equation 11) of   for HT ISR , given by The vectors q and v  are composed of the sensitivities of the QoI and observations, respectively, and capture all possible dynamical mechanisms via which perturbations in the control variables can change the QoI and observations. Therefore, DPP  (Equation 15) evaluates the degree to which HT ISR and   share their adjustment physics. Since wind adjustments dominate HT ISR and (most) observations   (Figures 4f-4j), a comparison of the τ y sensitivity map in Figure 4a with the τ y sensitivity maps in Figures 4b-4e elucidates the dynamical mechanisms and teleconnections that generate proxy potential. These mechanisms are discussed in detail in Loose et al. (2020) and only briefly reviewed in the next paragraph.
In Figures 4a-4e, sensitivities emerge in two main regions: (I) in the coastal wave guide along the eastern boundary of the subtropical North Atlantic; and (II) in topographic wave guides in the northeast Atlantic and the Nordic Seas (see Figures 5a-5c for labeled regions). Wind forcing in region (I) drives a pressure adjustment mechanism (Jones et al., 2018;Loose et al., 2020) which alters both the ISR geostrophic transport and the Irminger Current (Figure 3), causing anomalies in HT ISR and in all temperature observations (Figures 4a-4e). Wind forcing in region (II) spins up an anomalous barotropic circulation around Iceland (Loose et al., 2020), which simultaneously alters the Norwegian Atlantic and East Greenland Currents (Figure 3), causing anomalies in HT ISR and in the temperature observations θ A , θ C , θ D (Figures 4a, 4b, 4d and 4e). HT ISR and subsurface temperature in Denmark Strait, θ D , are the quantities that are most sensitive to the latter mechanism, driven by wind forcing in region (II). Consequently, θ D has the largest DPP for HT ISR among all observations considered (30%). The lower DPP of θ A (19%) is explained by the strong sensitivity of θ A to wind forcing in both regions (I) and (II). This results in destructive interference of information propagation because wind forcing in region (I) causes responses in HT ISR and θ A of equal sign, while wind forcing in region (II) causes responses in HT ISR and θ A of opposite sign (Figures 4a and 4b). The DPP of θ B and θ C is only 1%, since θ B is not sensitive to wind forcing in region (II) (Figure 5c),   (Table 1) to the magnitude of the vectors (f) q, . The yellow labels show the dynamical proxy potential (DPP) of   for HT ISR (Equation 15).

Obs
Location is mostly sensitivity to local forcing (Figure 5d). These characteristics lead to small sensitivity overlaps of HT ISR with θ B and θ C .

Data Redundancy versus Complementarity
Next, we demonstrate how DPP generalizes when considering multiple but still noise-free observations. We begin by combining θ A with either of the remaining temperature observations   , , , The observations θ A and θ B contain independent information. Indeed, θ A is sensitive to wind forcing in both regions (I) and (II) (Figure 5a), whereas θ B is sensitive to wind forcing only in region (I) (Figure 5b). Viewed in the {θ A , θ B }-informed subspace of the control space (Figure 5g), independent sensitivity information corresponds to the vectors v A and v B being close to orthogonal, with an enclosed angle of β = 74°. Consequently, is only a slight modification of v B (Figures 5b, 5d and 5g). In contrast to θ B , the observation θ D shows very similar sensitivity to wind forcing as θ A , concentrated in both regions (I) and (II) (Figures 5a and 5c). Similar sensitivity information is reflected by an angle of only δ = 30° between v A and v D (Figure 5h). Orthonormalization of the vector pair results in a removal of sensitivity to wind forcing in regions (I) and (II) (Figure 5e). The independent sensitivity information extracted from the observing systems ( Figure 5) is independent of the chosen QoI.
We now return to our QoI, and use the independent sensitivity patterns obtained in Figure 5 to quantify the DPP of the system  

{ , }
A  for HT ISR . DPP is obtained through a projection of q onto the  
The first term on the right hand side is equal to 19% (cf. Figure 4b)  The consequence of orthonormalization now becomes apparent. When combining θ A and θ B , data complementarity results in a DPP of (19 + 6)% = 25% (Figure 6b), exceeding the sum of DPP A = 19% and DPP B = 1% (Figure 6g). Put differently, considering θ A and θ B in combination helps to extract some of the observations' sensitivity information which is lost in destructive interference when treating θ A or θ B in isolation. The phenomenon of destructive interference was discussed in Section 3.2.1 and in Loose et al. (2020). Data complementarity can be viewed from yet another angle, when inspecting the position of q projected onto the {θ A , θ B }-informed subspace, denoted by q obs (Figure 6j). q obs is not aligned with either of the vectors LOOSE AND HEIMBACH 10.1029/2020MS002386 13 of 25 v A or v B . Therefore, a true information gain in the QoI is achieved when combining the observations θ A and θ B : the length of q obs increases, when advancing from the θ A -informed and θ B -informed components (orange dots, Figure 6j) to the {θ A , θ B }-informed component (radius of orange circle, Figure 6j).
Adding the observations θ D or θ C to θ A involves a certain degree of data redundancy, which is quantified in Figures 6h and 6i. Proxy potential of θ A for HT ISR originates in wind forcing in regions (I),(II) (Figure 5a); this sensitivity information is already contained in θ D (Figures 5c and 6k). Consequently, DPP(θ A , θ D ; HT ISR ) does not exceed DPP D = 30% (Figures 6c and 6h). Similarly, the relevant sensitivity information contained in the Irminger Sea surface observation θ C is already contained in the Irminger Sea subsurface observation θ A (Figure 6l). Thus, θ C does not lead to a gain in DPP when added to θ A (Figure 6d).
Finally, we are interested in the maximum achievable DPP for HT ISR , obtained by combining all four observations in our case study. Viewed within the three-dimensional subspace that is informed by the observing system {θ A , θ B , θ D } (Figure 5f), the {θ A , θ D }-informed yellow plane is almost orthogonal to the {θ A , θ B }-informed green plane (where the black plane is exactly orthogonal to the green plane). Hence, when adding θ D to the observing system {θ A , θ B }, the gain in DPP (green-yellow hatched, Figure 6e) almost coincides with Figure 6c), leading to a total DPP of 35% (Figure 6e). Completing the observing system by θ C does not increase the DPP any further (Figure 6f).

Noisy Observations
So far, our analysis has assumed noise-free observations. Next, we study the EPP of our observations   ; this notion does account for observatio  nal noise, as encoded in the noise covariance matrix R (Equation 1). Recall that the EPP of   for HT ISR is equal to the relative uncertainty reduction in HT ISR that is achieved when adding   to the underlying state estimation framework (Equation 12). Following the common assumption of uncorrelated observation errors (e.g., Forget et al., 2015), we only need to specify the diagonal entries of R, that is, the error variance  2  of each observation   . We assign   = 0.1°C for all observations (Table 2). We also consider the impact of varying ɛ C , by testing for ɛ C = 0.2°C and ɛ C = 0.3°C. The rationale for this addition is that climatological surface temperature, measured by θ C , is more variable than climatological subsurface temperature (Locarnini et al., 2013), and can therefore be expected to be more noisy.

Sensitivity-To-Noise Ratio
The strength of the constraint provided by each individual observation   is quantified by the eigenvalue   (Equation 5) corresponding to the   -informed direction v  (Figure 2a). It is given by  to atmospheric forcing than subsurface temperature (  Table 2), the SensNR of θ C is higher than that of θ A , θ B , θ D (Figure 7). This remains true if the noise variance for θ C (i.e.,  2 C ) is assumed four-or even nine-times as large as that of the subsurface observations (Figure 7).
is equal to the prior uncertainty in the observed quantity   (cf. Equation 8), i.e. the uncertainty given the prior knowledge in the ocean state estimate, before taking the actual measurement. Thus, observations with SensNR smaller than 1 (here: θ A , θ B , θ D , gray rectangle in Figure 7) are characterized by a prior uncertainty,  2 ( ) B  , that is smaller than their assumed observational uncertainty,  2  .
The EPP of   for HT ISR is given by  The factor   indicates what fraction of DPP  can be retrieved and will therefore be referred to as the "effectiveness" of the observation   . Note that, in contrast to DPP  , the observation's effectiveness,   , is independent of the QoI under consideration. Instead, it is solely determined by the observation's SensNR (Figure 7), observations with higher SensNRs are more effective. Therefore, the effectiveness of the surface observation θ C is higher than that of the subsurface observations θ A , θ B , θ D (Figure 7). In fact, the effectiveness of θ A , θ B , θ D is less than 50%, due to their SensNR being less than 1 (gray rectangle, Figure 7).
In this section, we studied the SensNR,   , and associated effectiveness,   , of each individual observation   . In the next section, we will establish a connection between the   and the set of eigenvalues   4 1 { } i i , where the latter set characterizes the observing system that is jointly formed by {θ A , θ B , θ C , θ D }.

Combining Noisy Observations
We now combine all four temperature observations of our case study, while taking into account their observational noise. In Equation 1, the resulting observing system is represented by LOOSE AND HEIMBACH 10.1029/2020MS002386 16 of 25  Table 2. An observation that falls into the light gray rectangle has sensitivity-to-noise ratio smaller than 1.
where the latter denotes a diagonal 4 × 4 matrix, with diagonal entries equal to the noise variances  2  , chosen as in Table 2. For the sake of brevity, we focus on the case ɛ C = 0.2°C. (Cases with alternative choices for ɛ C are presented in Figure S1). We compute the eigenvectors and eigenvalues,   By definition, the first eigenvector v 1 points in the direction of maximal curvature of  misfit ( ) J u . This direction is almost aligned with the θ C -informed direction, spanned by v C (Figure 8a), because the surface observation θ C has a much higher SensNR than the remaining observations (Figure 7). The remaining eigenvectors, v 2 , v 3 , v 4 , have little contribution from θ C (purple dots, Figure 8a), and are instead a linear combination of v  ,  , , A B D  (Figure 8b). The τ y component of v 2 (Figure 8d) extracts the dominant sensitivity patterns along the eastern boundary of the North Atlantic (region I), shared by θ A , θ B , and θ D , and in the northeast Atlantic and the Nordic Seas (region II), shared by θ A and θ D (Figures 5a-5c). The τ y component of v 3 (Figure 8e  governed by sensitivities in region (I), as set by θ B . Lastly, the τ y component of v 4 (Figure 8f) is dominated by sensitivity dipoles local to the three observing sites (yellow dots), which emerge due to the effect of Ekman pumping.
The eigenvalues λ i (inserted in Figures 8c-8f) are closely linked to the SensNRs   (Table 2) of all observations   involved, via the following relations (Bunch et al., 1978): Each eigenvalue λ i determines the effectiveness, η i = λ i /(λ i + 1), of the eigenvector v i (insets in Figures 8c-8f and diamonds in Figure 9b). We consider how the effectiveness of the observing system components changes with observational noise, and inflate the observational noise covariance by a factor α, to αR. As α varies from 0 (no noise) to 1 (full noise), effectiveness decays as λ/(λ + α), from 100% to λ/(λ + 1) (Equation B1).
Here, λ is a placeholder for either a SensNR   (Figure 9a) or an eigenvalue λ i (Figure 9b)  effectiveness of the surface observation θ C (Figure 9a) as well as v 1 (Figure 9b) is slower than that of the remaining observations and eigenvectors.
Decay in effectiveness causes decay in proxy potential. When considering individual observations   , inflating observational noise leads to a decay in proxy potential according to  (Figure 9c).
If instead, the full observing system is considered jointly, the decay in proxy potential is given by see Figure 9d (and Appendix B). The expression in Equation 21 involves the eigenvectors and eigenvalues of the misfit Hessian. For noise-free observations, proxy potential is equal to DPP (α = 0, pentagons in Figures 9c and 9d, Equation 11). It decays to EPP for fully inflated noise (α = 1, squares in Figures 9c  and 9d, Equation 12). Even though the surface observation θ C has highest SensNR and, thus, slowest decay in effectiveness (Figure 9a), its proxy potential for HT ISR is lower than that of the subsurface observations θ A and θ D , due to its almost negligible DPP at the very outset α = 0 (Figure 9c). Through a similar argument, v 1 contributes less to proxy potential than v 2 and v 3 (Figure 9d), despite its relatively highest effectiveness ( Figure 9b). The insignificance of θ C implies that proxy potential of the observing system {θ A , θ B , θ C , θ D } for HT ISR is essentially insensitive to the choice of the observation error ɛ C (Figures S1g-S1i).

Discussion
Hessian UQ and optimal observing system design have remained underexplored computational tools in oceanography, despite their successful application by the NWP and CSE communities. In this paper, we provided dynamical insight into Hessian UQ and how to leverage this method to design ocean climate observing systems. Our results warrant some general, conceptual remarks (Section 4.1), as well as discussion of specific implications for our North Atlantic case study (Section 4.2). We conclude with a discussion of dimension reductions of the Hessian (Section 4.3), limitations (Section 4.4), and an outlook (Section 4.5).

Conceptual Considerations
In the context of Hessian UQ, optimality refers to a well-defined goal of the observing system, often expressed in terms of one or several QoIs to be monitored (e.g., ocean transports) or predicted (e.g., sea-ice area). Given such a QoI, we rephrased the degree of optimality of an observing system in terms of 'proxy potential' (Equation 12), defined as the reduction in QoI uncertainty (on a scale of 0%-100%) that would be achieved if the observing system was added to the ocean state estimate. We showed that proxy potential combines three main concepts: (i) the degree of shared adjustment physics between QoI and observations, measured by the projections q•v i (Figures 4 and 8); (ii) data redundancy versus complementarity of the distinct members of the observing system, through orthogonality of {v i } (Figures 5 and 6); and (iii) the sensitivity-to-noise ratios (SensNR) of the observational assets, which determine the effectiveness of the observing system, through multiplication by the factors λ i /(λ i + 1) (Figure 9).
Concept (i) can be interpreted as the dynamical analog of statistical correlation between QoI and observations, where proxy origins (in the space of uncertain control variables) are unambiguously identified by the adjoint model (Loose et al., 2020).
Concept (ii) is associated with the eigenvectors v i of the misfit Hessian. Given their orthogonal structure, the eigenvectors may be compared to statistical empirical orthogonal functions (EOFs), but with the following important distinction. Whereas EOF patterns are based purely on statistical evidence, the Hessian eigenvector patterns are computed through the model's linearized dynamics. The eigenvector patterns arise as a linear combination from sensitivities of the observed quantities that are part of the observing system under investigation. They elucidate the dynamical connections between changes registered by the observing system and (local and remote) oceanic perturbations back in time that cause these observed changes ( Figure 8).
LOOSE AND HEIMBACH 10.1029/2020MS002386 19 of 25 Concept (iii) is linked to the eigenvalues λ i of the misfit Hessian. The eigenvalues quantify the strength of the constraints imposed by the observing system, and are determined by the SensNR of all observations that are part of the observing system. The SensNR of an observation can be regarded as the dynamical analog of the statistical signal-to-noise ratio (SNR). The SensNR accounts for all possible forcing scenarios (as captured by changes in the control variables), while the SNR is based on the statistics of observed or simulated data, which samples only certain forcing occurrences. The leading eigenvectors inherit the sensitivity patterns from observations with highest SensNR (Figures 8 and S1). The higher the observations' SensNRs (or eigenvalues), the more efficiently the observations reduce uncertainty in the state estimate and QoIs (Figure 9).

Implications for Our Case Study
We found that wind anomalies in the coastal and topographic wave guides, along the eastern and northern boundaries of the North Atlantic, are the largest origins of proxy potential (Figure 4). Such wind anomalies drive pressure adjustment mechanisms with a basin-wide response in North Atlantic circulation and temperature (Loose et al., 2020): adjustments that are registered by QoIs and observations alike, even if the observing system is remote from the QoI. This result can be rephrased in terms of UQ ( Figure 1) as follows. North Atlantic temperature observations reduce uncertainties in various atmospheric forcing variables at various locations (via path UQ1). Among these, it is primarily uncertainty reduction in surface momentum fluxes along the eastern and northern boundaries of the North Atlantic which leads to uncertainty reduction in the QoI, heat transport across the Iceland-Scotland ridge (via path UQ2).
For a given subsurface temperature observation in the Irminger Sea, our Hessian UQ analysis reveals that a subsurface temperature observation off the Portuguese coast would provide more independent information than a subsurface temperature observation in Denmark Strait (Figures 5 and 6). This result is explained by the fact that the Irminger Sea and Denmark Strait subsurface observations have very similar adjustment physics. Therefore, their information content is to a certain degree redundant.
Finally, our case study suggests that surface temperature observations have lower proxy potential for remote QoIs than subsurface temperature observations, despite their higher SensNR (Figure 9). This is due to strong sensitivity of surface temperature to local air-sea fluxes, which overrides their sensitivity to the large-scale, basin-wide adjustment mechanisms that are relevant for remote QoIs (Figure 4).

Dimension Reduction
One of the main computational challenges to UQ is the curse of dimensionality (Bellman, 1957). Since the uncertain parametric model inputs (or control variables) are typically adjusted on a grid point basis of the GCM, their number is large: . The calculation of the full Hessian-a matrix with 12 (10 )  to 16 (10 )  elements-would require years of extensive computer resources, an intractable endeavor.
The success of Hessian UQ relies on approaches that are more computationally efficient, two of which we consider: first, an a-priori reduction, and second, a data-informed reduction of the control space dimension. In this paper we have pursued the second approach, as we will further discuss in the next paragraph. In contrast, Kaminski et al. (2015Kaminski et al. ( , 2018 follow the first approach, by aggregating and adjusting their control variables uniformly over large regions (e.g., Figure 2 in Kaminski et al., 2015), rather than on a model grid point basis. This "large region approach" reduces their control space to a total of about 150 control variables, and it is then feasible to explicitly compute the full Hessian (150 2 entries). In practice, the large region approach requires to spatially accumulate sensitivities of QoIs and observed quantities over the pre-defined large regions, as exemplified in Figure 10. The eight regions defined in Figure 10c would reduce the dimension of our control space from 6 (10 )  (Table 1) to 8 ⋅ 4 = 32. However, the spatial accumulation of sensitivities would imply that proxy origins and adjustment mechanisms, for example, along the basin boundaries, are no longer resolved (Figures 10d  and 10e). Proxy potential would be artificially lost (right yellow label, Figure 10). Note that for other LOOSE AND HEIMBACH 10.1029/2020MS002386 20 of 25 QoIs, this approach could overestimate (rather than underestimate) proxy potential and uncertainty reduction.
Because of the ad-hoc nature of a-priori control space reductions, and the difficulties it incurs (Figure 10), we advocate the approach of data-informed reductions of the control space for the following reason. Even though the Hessian in our North Atlantic case study consists of 12 (10 )  entries (Section 3.1), the misfit Hessian is only of rank 4, equal to the number of observations involved. The four Hessian misfit eigenvectors with nonzero eigenvalues capture the Hessian's full information. They were extracted efficiently while preserving the physical mechanism that led to uncertainty reduction. The concept of data-informed control space reduction generalizes to large, complex observing systems, for example, mixed mooring arrays and autonomous instruments, which include thousands to millions of observations in time and space. While it becomes intractable to compute all (thousands to millions of) misfit Hessian eigenvectors, randomized numerical linear algebra for low-rank approximations can be used to extract the leading eigenvectors with highest eigenvalues (M′ ≪ M in Equation 3, Bui-Thanh et al., 2012;Kalmikov & Heimbach, 2014;Liberty et al., 2007). Moore et al. (2017) used a related technique in a regional ocean setting. They derived data-informed reduced-rank approximations of the Hessian, but with reductions sought in the observation space, rather than the control space. The two approaches are equivalent (or "dual" to each other), and the implementation of the underlying variational data assimilation scheme may determine which of the two approaches is more convenient to employ. We argue that an eigen-decomposition in the control space, as suggested here, has the appeal of a straightforward dynamical attribution of proxy origins.

Limitations
Some shortcomings of the method presented should be acknowledged. Hessian UQ relies on an accurate specification of the prior and noise covariance matrices, B and R (Equation 1). This is emphasized, for instance, by the fact that the relative weight of surface versus subsurface observational noise determines the observations' relative effectiveness, and thus the patterns that dominate the leading eigenvectors of the misfit Hessian ( Figure S1). A second shortcoming is that the results may suffer from model dependency, a problem common to all methods for model-informed observing system design. A third limitation is that Hessian UQ makes a Gaussian approximation of the posterior probability distribution for the uncertain control space and the estimated ocean state space. This approximation is accurate if the linearized model   (Table 1) to only 32, via the "large region approach." The yellow labels show the values for dynamical proxy potential (DPP; Equation 11), using the full control space (left, cf. Figure 4b) and the a-priori reduced control space (right). On the right, DPP evaluates to 0% due to an entire cancellation in the projection q • v A . For instance, sensitivities of equal sign in the subtropical Atlantic (region 2, (d) versus (e)) make a strongly positive contribution to q • v A , while sensitivities of opposite sign in the Norwegian Sea (region 5, (d) versus (e)) make a strongly negative contribution.
provides a good representation of the ocean dynamics on the time scales investigated. The results by Loose et al. (2020) indicate that on the five-year time scale considered, nonlinearity is not a major obstacle, at least not in the non-eddy-resolving model under consideration. In situations where strong nonlinearities are barriers to Gaussian approximations, Hessian UQ in combination with non-Gaussian sampling methods have shown promise (Petra et al., 2014), but are yet to be explored in ocean and climate modeling. For instance, Stochastic Newton MCMC employs a local Gaussian approximation (given by the local inverse Hessian), which is then used as a proposal distribution for the posterior probability distribution (Petra et al., 2014).

Outlook
In our case study, we made simplifying assumptions regarding the control variables and the prior error covariance matrix (Table 1) to enable a clearer understanding of the methodology. These simplifications are readily relaxed in future work. Based on the insights gained here, we aim to compute reduced-rank approximations of the Hessian for large observing systems within the ECCO framework. Our case study highlights that the stopping criterion for truncating the eigenvalue spectrum has to be chosen carefully, because the leading Hessian eigenvectors are not always the most important ones for informing a given QoI. Indeed, eigenvectors lower down in the spectrum capture important dynamical teleconnections originating from the sensitivity of subsurface (rather than surface) observations. Future work should address the interesting question whether the abundance of surface observations (available from satellite altimetry) and their mutual complementarity (due to their local sensitivity) may be able to cover for the large-scale sensitivities of subsurface observations. The technique presented in this paper is complementary to the more widely used OSSEs. Hessian UQ elucidates dynamical teleconnections that communicate observational constraints-via ocean currents, wave dynamics, Ekman dynamics, and geostrophy-over basin-scale distances and on monthly to interannual time scales. It provides an approach for guiding the design of observing systems that (i) maximize the information about (possibly remote) QoIs that are difficult or impossible to observe directly, and (ii) are complementary to the existing observational database. We hope that Hessian UQ, in combination with OSSEs and other tools, will be more widely used for tackling the grand community challenge of co-designing a cost-effective and long-term Atlantic observing system in the coming years.