Progress in Oceanography The accuracy of estimates of the overturning circulation from basin-wide mooring arrays

Previous modeling and observational studies have established that it is possible to accurately monitor the Atlantic Meridional Overturning Circulation (AMOC) at 26.5°N using a coast-to-coast array of instrumented moorings supplemented by direct transport measurements in key boundary regions (the RAPID/MOCHA/WBTS Array). The main sources of observational and structural errors have been identi ﬁ ed in a variety of individual studies. Here a uni ﬁ ed framework for identifying and quantifying structural errors associated with the RAPID array-based AMOC estimates is established using a high-resolution (eddy resolving at low-mid latitudes, eddy permitting elsewhere) ocean general circulation model, which simulates the ocean state between 1978 and 2010. We de ﬁ ne a virtual RAPID array in the model in close analogy to the real RAPID array and compare the AMOC estimate from the virtual array with the true model AMOC. The model analysis suggests that the RAPID method underestimates the mean AMOC by ∼ 1.5Sv (1Sv=10 6 m 3 s − 1 ) at ∼ 900m depth, however it captures the variability to high accuracy. We examine three major contributions to the streamfunction bias: (i) due to the assumption of a single ﬁ xed reference level for calculation of geostrophic transports, (ii) due to regions not sampled by the array and (iii) due to ageostrophic transport. A key element in (i) and (iii) is use of the model sea surface height to establish the true (or absolute) geostrophic transport. In the upper 2000m, we ﬁ nd that the reference level bias is strongest and most variable in time, whereas the bias due to unsampled regions is largest below 3000m. The ageostrophic transport is signi ﬁ cant in the upper 1000m but shows very little variability. The results establish, for the ﬁ rst time, the uncertainty of the AMOC estimate due to the combined structural errors in the measurement design and suggest ways in which the error could be reduced. Our work has appli-cations to basin-wide circulation measurement arrays at other latitudes and in other basins as well as quantifying systematic errors in ocean model estimates of the AMOC at 26.5°N.


Introduction
Estimating ocean transports of volume, heat and freshwater is a fundamental oceanographic activity that provides data critical to the study of the ocean's role in the mean climate and climate variability. Bryden and Hall (1980) highlighted that the Atlantic Meridional Overturning Circulation (AMOC) is responsible for a large part of the climatically important ocean heat transport, which is the largest of any ocean basin (Bryden and Imawaki, 2001). The phenomenology, dynamics and climate impacts of the AMOC are a complex subject which despite extensive research are still to be fully understood (see Buckley and Marshall, 2016 for a comprehensive review of what is known and not known about the AMOC). The AMOC is typically defined as the maximum of the zonally (basin-wide) and vertically integrated meridional transport, which occurs near 30°N, coincident with the maximum poleward heat transport. The first estimates of the strength of the AMOC and associated heat transport were based on hydrographic sections (typically at 24.5°N e.g. Bryden et al., 2005) that provided a snapshot of the circulation but a growing awareness of the variability of the AMOC stimulated the development of continuous measurement using mooring arrays. The RAPID array (Cunningham et al., 2007) at 26.5°N in the Atlantic Ocean has provided continuous, basin-wide, fulldepth observational estimates of the AMOC since 2004.
The RAPID/MOCHA/WBTS array (hereinafter referred to as the RAPID array) has revolutionized basin scale oceanography by supplying continuous estimates of the meridional overturning transport , and the associated basin-wide transports of heat (Johns et al., 2011) and freshwater (McDonagh et al., 2015) at 10day temporal resolution. These estimates have been used in a wide variety of studies characterizing temporal variability of the North Atlantic Ocean, for instance establishing a decline in the AMOC between 2004 and 2013 (Smeed et al., 2014), recording of a substantial downturn in the AMOC in -10 (McCarthy et al., 2012 and its subsequent effects on North Atlantic heat content  and sea surface temperature  and winters in North Western Europe (Buchan et al., 2014, Maidens et al., 2013. The RAPID array has also been important in defining shorter timescale variability, including seasonal variations in the AMOC and their origin Chidichimo et al., 2010;Mielke et al., 2013;Duchez et al., 2014;Pérez-Hernández et al., 2015), and extreme volume and heat transports on sub seasonal timescales (Moat et al., 2016;Cunningham et al., 2013).
With such a large number of studies reliant on the RAPID measurements and other basin-wide array data, it has become imperative to accurately establish the structural error of the RAPID array estimates as has been done for the observational error (see McCarthy et al., 2015). The original proof of concept studies which established the feasibility of the RAPID array were conducted using the (at the time) high-resolution (1/4°-i.e. eddy permitting) OCCAM and FLAME (1/3°) ocean general circulation models (Hirschi et al., 2003Baehr et al., 2004;Hirschi and Marotzke, 2007). The demonstrated a correspondence between the true model transports and a proxy based on the basin-wide geostrophic transport, calculated from moorings on the eastern and western boundaries and over the Mid Atlantic Ridge (MAR). We do not seek to challenge these findings, but rather wish to put them on a firm quantitative footing for the real array, in all its complexity, using a state-of-the-art high resolution (eddy resolving) ocean general circulation model (OGCM).
A key problem when using geostrophic calculations is the need for an absolute velocity at the reference level (strictly, the true or absolute geostrophic velocity at the reference level). This issue has been recognized since the dawn of modern oceanography, often in the context of determining the volume transport through a hydrographic section and a number of methods have been developed to address it, from the assumption of one or more levels of no motion (e.g. Bryden et al., 2005), to much more sophisticated approaches, where further constraints such as a requirement for zero net volume transport are imposed (inverse methods, exemplified by e.g. Wunsch, 1978). Ganachaud (2003) provided a thorough exposition of potential error in one-time hydrographic sections when using the inverse method at 25°, 36°and 48°N, anticipating many of the concerns of this paper, in particular the reference level error was found to be of order ±3 Sv at 25°N. Unsampled regions ("bottom triangles") and ageostrophy were by contrast found to be much smaller, although measurement error due to internal wave processes were found to be of the same order of magnitude to the reference level error. Rather than using classical inverse methods per se, Hirschi et al. (2003) suggested the simple method of adding a spatially constant but time variable barotropic velocity to the measured basin-wide geostrophic velocity in order to ensure zero net volume transport across the basin. In an ocean with uniform depth this method of solution will always give the correct answer even if the reference velocity varies with longitude. However, when the water depth varies across the array, then the zero net-transport assumption is not in principle sufficient and variations in the vertical structure of the circulation solution dependent on the choice of reference level emerge . Whilst Hirschi et al. (2003) acknowledged this issue, they did not perform further analysis on their model simulations.
Several studies have highlighted the potential errors introduced by the imposition of the zero net-transport constraint (Searl et al., 2007;Hughes et al., 2013). Kanzow et al. (2009) and McCarthy et al. (2012) investigated the accuracy of this assumption by comparing the transport variability derived from basin-wide pressure differences in bottom pressure recorders with the transport variability derived from the application of the mass compensation constraint. The results of these studies suggest that the error associated with the mass constraint in the estimate of the maximum of the overturning stream function is comparable to the accuracy of the other elements of the array at 1.5 Sv . Further independent validation of the mass constraint was provided by Landerer et al. (2015) who verified the detrended interannual variation associated with the mass constraint using bottom pressure data derived from the GRACE gravity satellite.
There have been further concerns raised about the accuracy of the RAPID array, notably by Wunsch (2008) who worried that the dependence of the geostrophic transport calculation on the endpoints would result in large errors in the basin-wide transport estimate due to eddy variability on the boundaries. This argument has been countered by Kanzow et al. (2010) who showed that eddy variability rapidly decreases close to the boundaries, resulting in a much lower error than suggested by Wunsch's analysis.
The assumed low magnitude of the ageostrophic transport has been further questioned by Stepanov et al. (2016) who suggest that ageostrophic transport associated with mesoscale eddies in the western North Atlantic makes a significant contribution to the volume transport (and more particularly the heat transport). As the RAPID array is not able to capture this effect, it was suggested that RAPID may be missing a significant part of the meridional volume transport.
The operational design of the RAPID array introduced further assumptions, including concatenation of geographically separated moorings . These assumptions have never been rigorously scrutinized for the errors they have introduced into the volume transport estimates, although some authors have discussed the potential for errors (e.g. Roberts et al., 2013 in the eddy permitting model regime). Stepanov et al. (2016) compared model proxies (in the eddy resolving regime) calculated in a very similar way to that done using the real RAPID array, with model truth and found that the proxy slightly overestimates the model truth.
The various sources of error when using a RAPID-type array can be rigorously quantified in a model study as absolute geostrophic currents can be determined in the model given knowledge of model pressure. Ganachaud (2003) did in fact estimate the variability in the reference level velocity using absolute geostrophic currents from an eddy permitting ocean general circulation model. However he did not take the further step of conducting a perfect model study in order to directly compare an estimated geostrophic transport using the same method as with observations with the absolute geostrophic transport from the model. Furthermore, the RAPID array was not in place then, and therefore Ganachaud's work was not performed in the context of an operational array giving continuous estimates of basin-wide transport.
In this paper we go further and rigorously separate the structural bias into (sometimes mutually compensating) components due to reference level assumptions, unsampled regions, and ageostrophy.
Our analysis enables us to accurately pinpoint the magnitude, nature and spatial location of structural errors in the RAPID array-based AMOC estimates for the first time. This reaffirms the ability of the array to provide reliable transport estimates, and leads to new insights that can inform array design improvements and thus provide more accurate future estimates of the AMOC and heat transport (two key climate parameters).
In contrast to previous studies, we make use of the sea-surface height information provided by our free-surface model in order to determine the true or absolute geostrophic transport and thereby establish a rigorous framework to identify and quantify structural errors in the RAPID array and other arrays like it.
In this paper we examine sources of error due to: (i) use of a fixed reference level and the assumption of zero net B. Sinha et al. Progress in Oceanography 160 (2018) 101-123 transport (ii) neglect of ageostrophic terms (iii) neglect of regions not covered by the array instruments Using a high-resolution OGCM (1/12°-eddy resolving in low-mid latitudes, eddy permitting elsewhere), we first quantify the mean bias between the proxy method and the model truth and the time variable error on both sub-annual and inter-annual timescales in different depth ranges.
We develop a mathematical analysis which gives the error in the geostrophic transport in terms of the average velocity at the base of each mooring pair (similar to the analysis of Hughes et al., 2013) and test this in the OGCM. Additionally we use the model to estimate the magnitude of the transports in regions not covered by the array, and the ageostrophic transport and its spatial variation. Less significant as a source of bias, but nonetheless interesting is the assumption that the Ekman transport is distributed uniformly within an Ekman layer 100 m deep. We quantify the resulting bias in the transport estimates and verify that the classical formula for the Ekman transport accurately reflects the net Ekman transport, albeit not its vertical variation. Finally, we discuss the effect of assuming a constant value of the Coriolis parameter for each mooring pair, equal to the average of the individual values (since the latitude of deployment varies) and demonstrate that this makes a very small (ageostrophic) contribution to the bias.
We quantify the relative magnitude of these terms as a function of depth and also divide the bias into parts associated with each component of the monitoring array.
Lastly, we investigate the time variation of the bias terms and the relationship between biases at different depths. We discuss potential improvements in the RAPID array based on our results and also discuss which aspects of the results are likely to be robust given potential model biases.

Model description
We employ the NEMO (Nucleus for European Modelling of the Ocean) general circulation model (Madec, 2008) implemented on the tripolar ORCA12 horizontal grid with a nominal 1/12°horizontal resolution and 75 vertical levels ranging in thickness from 1 m at the surface to ∼200 m at depth. South of 20°N the ORCA12 grid has 1/12°r esolution in longitude and the latitudinal resolution adjusts to match. North of 20°N the grid distorts, but retains matching latitudinal and longitudinal resolution. At 26.5°N, the location of this study, the grid distortion is negligible and model grid lines coincide with lines of latitude, with a grid resolution of 8.6 km. The model uses the Louvain La Neuve (LIM2) sea ice model (Fichefet and Maqueda, 1997) and was initialized with World Ocean Atlas Locarnini et al., 2006) data. The model was forced with the Drakkar Forcing Dataset Set version 4.1 (version 5.1 from 2008 onwards) (Brodeau et al., 2010) which supplies surface winds, air temperature and humidity at 6 h intervals, longwave and shortwave radiative fluxes at daily intervals, and precipitation (rain and snow) at monthly intervals. These datasets are temporally and spatially interpolated to the model grid and time step. A weak relaxation of the surface salinity to climatological values was applied using a piston velocity of −33.333 mm/ day/psu equivalent to a ∼180 day restoration timescale over a typical 200 m ocean surface mixed layer. The simulation period was 1 January 1978 to 31 December 2010 (33 years). Further details concerning the model configuration can be found in Marzocchi et al. (2015) and in Appendix C.7, Table C.7.1.
The performance of the NEMO model (particularly in the North Atlantic and at the location of the RAPID Array) has been evaluated in a number of previous papers including those by Marzocchi et al. (2015) (North Atlantic subpolar gyre), Hirschi et al. (2013) (mesoscale variability), Deshayes et al. (2013) (freshwater transports), Moat et al. (2016) (heat transport), Duchez et al. (2014) and Blaker et al. (2015) (both for the AMOC). The model displays a good simulation of the subtropical and subpolar gyres, including a realistic Gulf Stream separation and North Atlantic current path. The model also simulates the interannual variability at 26.5°N well, with a good representation of the major dips in transports in 2009 and 2010 and their subsequent effects on subtropical heat content . We note that there is some evidence that higher resolution, eddy resolving models are better able to correctly reproduce the observed mean AMOC and associated meridional heat transport at 26.5°N . Of most relevance for the present study, Fig. 1b of the paper by Blaker et al. (2015) compares the AMOC as estimated by the RAPID array to the AMOC calculated using data from the same NEMO simulation as used in this paper, for the period 2004-2011. It shows similar comparisons for three components of the AMOC, namely the Florida Strait transport, the Ekman transport and the Upper Mid Ocean transport, the latter referring to the southwards transport between the Bahamas and the African coast, excluding the northwards Ekman transport in the very surface (up to a few tens of meters) layer. This shows that the model simulates the variability of the AMOC very well, but with a mean value of about 15 Sv it underestimates the mean observed AMOC value by a few Sv. Slightly higher AMOC values have been simulated in other high-resolution ocean models (Xu et al., 2014;Iovino et al., 2014;Stepanov et al. 2016). Whereas it is beyond the scope of the present paper to explore the precise reason for such differences is it worth pointing out that the models used in Xu et al. (2014), Iovino et al. (2014) and Stepanov et al. (2016) differ from our model set up in ways (e.g. choice of surface forcing dataset; use of strong restoring for temperature and salinity; type of vertical coordinate, sea ice model) that inevitably affect the AMOC solution the models will adopt. However, such differences are not a problem for our study. Our main goal is not to simulate a "perfect" AMOC (compared to observations) but to better understand sources of error inherent to the RAPID AMOC monitoring array. Despite there being a bias in the mean AMOC there is very good agreement between observed and simulated AMOC and of particular importance   B. Sinha et al. Progress in Oceanography 160 (2018) 101-123 for the present study, NEMO also simulates the individual components of the AMOC well. The Ekman transport is virtually identical between model and observations, as expected. The mean Florida Strait transport is a little low compared to observations and its variability is not well reproducedhowever we do not expect it to be given that the Florida Strait transport is affected by eddies and waves which will have different timings in the model compared to the real ocean. The Upper Mid Ocean transport is reasonably well captured in terms of mean and seasonal cycle, but intra seasonal and interannual variability shows some divergence from observations, probably due to the presence of mesoscale variability. Summary statistics, including mean values and correlation coefficients are presented in Table 2 of Blaker et al. (2015). However Blaker et al. (2015) did not compare the RAPID observations with model output chosen to be as closely analogous to the real moorings as possible.
We deploy virtual moorings in the model in close analogy to the RAPID array (see McCarthy et al. (2015) for details) and calculate the proxy transport in exactly the same way, as described below, with the following notable exceptions: (1) The mooring line is taken as a model grid line, which is zonal to a very good approximation, rather than the actual variable-latitude array used in the real world. The horizontal locations of the moorings are at the same longitudes as the real moorings, but displaced to the north or south in order to lie on the model gridline ( Fig. 1). (2) The reference level for geostrophic transports is taken as 4320 m compared with 4820 db (∼4750 m) for the real array ( Fig. 2a) because the boundary between North Atlantic Deep Water (NADW) and Antarctic Bottom Water (AABW) is higher up in the water column in the model compared with reality. (3) The western boundary mooring is taken as a single tall mooring at one location rather than as a concatenation of two moorings slightly displaced from one another (Fig. 2a). This is because the second mooring in the real array is to the east of the first and is entirely below 4320 m, the model reference level. Similarly there is one real mooring on the eastern boundary that lies entirely below 4320 m and this is not included in the virtual array. (4) Florida Strait, Western Boundary Wedge (WBW) and AABW transports (defined below, see also Fig. 2a) are taken as the true model transports and no attempt is made to sample these regions in a way analogous to the real world array (e.g. with virtual current meters at discrete locations in the WBW). (5) Because the real array records no data in the top 50 m (see Haines et al., 2013 introduced seasonal extrapolation to the surface using deeper measurements. No attempt is made to reproduce this in the model proxy calculations and all model output (18 vertical levels) in the top 50 m are used. (6) Model z-coordinate surfaces vary slightly as a function of sea-surface height, hence the model moorings sample the ocean at depths, which vary with time, unlike the real world moorings which sample at fixed depths (excepting knockdown by ocean currents). This feature of the model results in a small correction to the model pressure gradient in order to give the correct absolute geostrophic transport and for consistency, this correction is also added to the expression for the proxy transport. Further details are given in Appendix C.4. (7) A constant value for the Coriolis parameter is used, whereas the estimate from the real array uses a different value (depending upon latitude) for each mooring pair. Variations in Coriolis parameter along the section introduce a negligible ageostrophic transport: this is discussed further in Appendix A.

Theory
In this section we draw together the contributions of many earlier authors to build a comprehensive framework for diagnosing and attributing biases. Given the importance of the RAPID array it is important that a study specific to the RAPID array is conducted which exactly mimics the real arraysomething that has not been done thus far. Whilst focusing in detail on the RAPID array, there are clear applications to other basin-wide measuring arrays (e.g. Lozier et al. 2017;

RAPID proxy calculation
The objective of RAPID is to estimate the overturning streamfunction, ψ, where v(λ, z, t) is the meridional seawater velocity (zero by definition outside the ocean boundaries) along a line of constant latitude, φ, (positive northwards) as a function of longitude, λ, and depth, z (positive upwards). λ e and λ w are the eastern and western boundaries, a is the radius of the Earth and z′ is a dummy variable of integration.
The AMOC is defined to be the maximum of the overturning which occurs near 30°N in the subtropical North Atlantic at a depth of 900-1100 m. The RAPID array is nominally located along the φ = 26.5°N latitude line but in fact varies from 23.9°N to 27.9°N (red line, Fig. 1). In the NEMO ocean model we approximate ψ using model northward velocities (v) along a nearly zonal model grid line located at φ = 26.43°N and we define a model analogue to the RAPID observation-based proxy based on virtual moorings placed along the model gridline in a similar configuration to the actual RAPID array (black line, Fig. 1). We adopt a local zonal coordinate, x, representing the distance along the model grid line from x = X w at the western boundary to x = X e at the eastern boundary, v is reinterpreted as the northward velocity as a function of x, z and t, and ψ is reinterpreted as the streamfunction as a function of z and t. There are eight mooring locations analogous to the real RAPID array mooring locations: Figs. 2a,2b). Mooring 1 extends from the reference level Z REF to the surface, moorings 2, 3 and 4 extend from Z REF to just above the maximum height, Z MAR , of the Mid Atlantic Ridge (MAR) and moorings 5-8 between them cover the remainder of the water column up to the surface. The deepest depth of the ocean is Z -H . Thus Z 1 -Z 4 , are all taken as the reference level, Z REF , and Z 5 , coincidentally is equal to Z MAR . (Table C.1 lists the positions and depth ranges of the moorings and their correspondence with the real moorings). We define sea-surface heights at the mooring locations, The positions of the Florida Strait (FS), the Western Boundary Wedge (WBW) and seven uninstrumented (unsampled) regions, 1,2, ,7 i are demarcated as are seven rectangular regions where geostrophic velocities are calculated from pairs of moorings. In particular, for reasons that will become clear in Section 3.4, the geostrophic velocities at the bases of these regions , , , , , , 21 32 43 54 65 76 87 , are marked.
As comprehensively described by McCarthy et al. (2015), the RAPID method calculates the total transport per unit depth (i.e. the z-derivative of the streamfunction ψ), T(z,t), as the sum of (i) the Florida Strait transport (ii) the Western Boundary Wedge transport (both directly measured), (iii) the wind driven Ekman transport (estimated from the surface zonal wind stress), (iv) the geostrophic transport between X 1 , and X 8 with respect to a deep reference level, (v) an assumed time constant Antarctic Bottom Water flow below the deep reference level and (vi) a compensation term to satisfy a mass constraint of no net flow across the section. Model equivalents of these are as follows: Florida Strait transport, T FS , is calculated from model prognostic velocities in Florida Strait (hence neglecting uncertainties arising from the cable measurement calibration).
where X b is the location of the Bahamas (Fig. 2a). Recall that sea water velocity is defined to be zero in non-oceanic regions. Western Boundary Wedge transport, T WBW, is also calculated from model prognostic velocities (neglecting uncertainties related to sampling by the current meters and in particular extrapolation to solid boundaries) Ekman transport, T EK, is calculated from the model zonal wind stress, τ x : where A z ( ) is the width of the ocean basin excluding the Florida Straits and The Ekman velocity is assumed to be uniform in an Ekman layer of constant depth Z EK = −100 m. The Coriolis parameter, ∼ f , is assumed constant and equal to the average value along the model gridline on which the virtual moorings are deployed and ρ 0 is a reference density.
For the geostrophic transport, our model version of the RAPID method uses in situ density from the tall mooring at the western boundary, the moorings either side of the Mid Atlantic Ridge and an "eastern boundary mooring" constructed by vertically concatenating in situ density profiles from the five short moorings located along the eastern boundary continental slope. Thus the method uses four "moorings" at positions X 1 , X 2 , X 3 , and X 4 , which we refer to as M 1 , M 2 , M 3 , and M 4 (Fig. 2a). The dynamic height at depth z for the m th mooring, D m , is calculated by integrating in situ density for that mooring, ρ Mm from a reference level depth of Z REF = −4320 m to depth z: where g is the acceleration due to gravity. It is important to understand the difference between ρ 1 , … ρ ρ , , 2 8 which are the in situ densities at the actual mooring sites, and ρ M1 , M 2 4 which are used for the proxy transport calculation. ρ M1 , ρ M2 , ρ M3 are identical with ρ 1 , ρ 2 , ρ 3 but ρ M 4 equals ρ 4 below z = Z 5 = Z MAR ; ρ 5 between Z 5 and Z 6 ; ρ 6 between Z 6 and Z 7 ; ρ 7 between Z 7 and Z 8 ; and ρ 7 between Z 8 and the surface. The estimate for the geostrophic transport per unit depth with respect to a fixed reference level (before mass compensation) is defined as Note that the presence of the Mid Atlantic Ridge splits the Atlantic basin into two parts below a depth level, Z MAR . At shallower depths a correction is applied to the dynamic height difference across the basin to ensure continuity of dynamic height at depth Z MAR (the third and fourth terms on the RHS of the top line of Eq. (7)). A small extra correction term due to the variation of model coordinate surfaces in time is added to the RHS of Eq. (7) as described in Appendix C.4.
The AABW transport, T AABW , is calculated from model prognostic velocities below the deep reference level: where the overbar denotes a long term time average over year 1990-2004 of the model simulation (avoiding the model spinup period). As in RAPID, this specifies a time-invariant transport per unit depth below the reference level. The compensation term, C(z,t) = c(t)A(z) is chosen such that a basin-wide constraint of zero net transport is satisfied as for the real RAPID array (c(t) is a constant compensation velocity and A(z) is the basin width introduced earlier): is the total area found by integrating the width, A z ( ), of the depthlongitude section excluding the Florida Strait, but including the WBW. The compensation term is dependent upon all of the other terms.
Putting all the terms together (Eqs. (2)-(4) and (7)- (9)) we obtain our estimate of the stream function ψ p : Deeper than depth Z REF only the compensation and AABW terms contribute to ψ p . Appendix C gives details of how the various transports described here are calculated using model diagnostics.

Decomposition of the true model streamfunction
In order to evaluate the bias in the estimate, we decompose the actual streamfunction ψ as follows: The absolute geostrophic transport, T G , is obtained by integrating the local pressure gradient divided by the section-average Coriolis parameter from the western mooring to the eastern mooring at each depth. Applied between each of the moorings as illustrated in Figs. 1 and 2, taking into account the presence of the MAR, allows us to express the absolute geostrophic transport, T G , across the section as: is the horizontal density difference between moorings m and n, and is the sea-surface height (SSH) difference between moorings m and n.
. A small extra correction term due to the variation of model coordinate surfaces in time is added to the RHS of Eq. (13) as described in Appendix C.4. The ageostrophic transport, T AG , is given by Thus T AG includes the (true) Ekman compensation as well as any other ageostrophic flow (apart from that due to the variation of the Coriolis parameter along the line of the array). T Δ is the total of all the transport in all the unsampled regions (above Z REF ), described fully in Section 3.5.

Decomposition of bias in proxy estimate
Subtracting the full streamfunction (Eq. (12)) from the RAPID-style estimate (Eq. (11)) gives We have thus identified four terms that contribute to the bias in the proxy. These are related to (i) the mismatch between the actual and proxy geostrophic transports, due to an unknown reference level (T P -T G ) (ii) ageostrophic flow (T EK -T AG ) (iii) flow in unsampled regions T AABW − T AABW − T Δ and (iv) the mass compensation term C. We can split the mass compensation term in order to individually compensate each of the other three terms: The residual in Eq. (19) arises because the depth integral of the proxy transport is set to zero, Eq. (9), hence the depth integral of the bias is non-zero and part of the compensation must cancel the depth integrated bias.

Further decomposition of biases due to reference level error
It can be shown (see Appendix B) that the error due to the unknown reference level for the calculation of geostrophic transport can be related to the geostrophic velocities at the bases of the mooring pairs. At any given depth level, all mooring pairs below that depth will contribute to the bias and conversely any mooring pairs above that depth will not contribute. The contribution is simply the geostrophic velocity at the base of the mooring pair multiplied by the horizontal distance between the moorings and the vertical distance across the base of the mooring pair: where where ρ Δ mn and h Δ mn are as defined previously (Eqs. (14) and (15)

Further decomposition of biases due to unsampled regions
The total transport per unit depth in the unsampled regions is the sum of the seven bottom triangles illustrated in Figs. 2a and 2b: where the transport in each bottom triangle is estimated by . For technical and numerical details of how these biases were calculated from the model output see Appendix C.

Mean bias
We first examine the estimated streamfunction, ψ p , in comparison to the model truth ψ (Fig. 3a). The true streamfunction, calculated over 15 years (years 1990-2004) to avoid the model spinup period, has a maximum of 15.2 Sv at 26.5°N at a depth of ∼857 m. The net flow through the section is 1.6 Sv southwards, composed of the Bering Strait outflow in the model plus any net loss or gain of volume due to evaporation and precipitation between 26.5°N and the Bering Strait. The mean proxy streamfunction displays a realistic shape, including a maximum value at the correct depth, however it underestimates the true streamfunction at virtually all depths, in particular at the depth of the maximum, ∼857 m the mean value of the proxy streamfunction is 13.6 Sv, so there is a mean bias of −1.6 Sv. The proxy estimate of the net flow was set to zero and so is different when compared to the true net flow (−1.6 Sv) over years 1990-2004. This results in a proxy estimate of the AABW transport (below 4320 m) that is slightly lower than the true streamfunction. Additionally, the proxy AABW transport is reduced by the mass compensation. It would be possible to require that the net throughflow of the proxy estimate, instead of being zero, be set to the actual throughflow of −1.6 Sv. At ∼ 1000 m depth the effect of this would be to increase the bias by a small amount (∼0.1 Sv, not shown).
Also shown is the streamfunction based on the absolute geostrophic transport (T ) G obtained from the model pressure gradient, not including the unsampled regions (see Appendix C.4), which has a maximum value of 14.7 Sv at 857 m, in between the other two curves. The streamfunction based on absolute geostrophic velocities need not show zero net transport because there may be balancing transport in the ageostrophic flow or the un-sampled regions.
The difference between the geostrophic streamfunction and the proxy streamfunction (red and black lines respectively in Fig. 3a) therefore reflects only the reference level bias and the mass compensation. On the other hand, the difference between the geostrophic and the true streamfunction (red and green lines in Fig. 3a) reflects only the unsampled regions and the ageostrophic transport. The geostrophic streamfunction corresponds closely with the true streamfunction down to about 3000 m depth which suggests that the unsampled regions and ageostrophic currents are relatively small in this depth range (although it is possible that they cancel each other) and most of the difference of the true streamfunction with the proxy streamfunction must be due to the reference level term and its compensation. Between 3000 m and 4320 m there is a much larger difference between the true streamfunction and the geostrophic streamfunction and this implies a large contribution from terms other than the reference level bias.
The relationship between the proxy, geostrophic and true streamfunction is further illustrated in Fig. 3b by removing the Ekman transport (Eq. (4)) from each quantity, focusing only on the instrumented regions (i.e. removing T Δ from the true streamfunction). The green and red curves now overlie to a large extent even at depths deeper than 3500 m and in particular show a very similar net throughflow of ∼−3Sv. The difference between these curves is purely due to ageostrophic transport. For comparison with Fig. 3a, the maximum values at 857 m depth are 12.5 Sv for the true streamfunction minus Ekman and unsampled regions, 10.4 Sv for the proxy minus Ekman, and 11.5 Sv for the absolute geostrophic transport minus Ekman. The time series of the streamfunction minus the Ekman component are shown in Fig. 3d. This clearly shows that time variations of the bias are relatively small (on annual timescales) and do not compromise the ability of the method to resolve variability of the AMOC on interannual and decadal timescales. The relative contributions of the various terms are calculated more precisely in Section 4.4.

Ekman transport
In reality, and in the model, the Ekman transport has a depth dependence. This is illustrated by Fig. 3c, which shows the total transport per unit depth and its decomposition into geostrophic and ageostrophic components in the top 120 m of the water column. Ekman theory (which assumes a constant vertical viscosity) suggests that the Ekman velocity should be at a maximum at the surface and weaken exponentially with depth. In the model, the northwards Ekman transport per unit depth (approximately identical in the surface layer with the ageostrophic velocity) has a peak of ∼0.073 Sv/m at ∼5 m depth and reduces to small values at 100 m depth (black line in Fig. 3c). There is a good correspondence in the model between the net ageostrophic transport in the top 100 m (total-geostrophic) and the Ekman transport estimated from the surface wind stress (respectively 3.4 Sv versus 3.2 Sv averaged over years 1990-2004, not shown). By 100 m the ageostrophic velocity remains roughly constant at ∼0.001 Sv/m.

Accuracy of the proxy streamfunction at different depths and timescales
In order to assess the accuracy of the RAPID proxy at different timescales and depths, we present scatter plots of proxy versus true streamfunction (both detrended) for unfiltered, and interannual timescales (Fig. 4). Considering the raw 5 day means output by the model (Fig. 4a) the proxy method achieves a squared correlation coefficient of 0.9 for 0-900 m volume transport, explaining 90% of the variance. The bias is approximately −1.5. Sv (relative to the 15.2 Sv model mean value) and the accuracy (RMS error) of the correlation is 1.2 Sv. At interannual timescales (Fig. 4d) the variance explained is similar and the accuracy is much higher (0.29 Sv RMS). On the other hand at mid depths (2022 m) for unfiltered data, the bias is ∼0.8 times the true streamfunction and the RMS error is 1.7 Sv. The correlation is weaker, though still high for mid depths (Fig. 4b, e) and very weak for the deepest depths (Fig. 4c, f). The reduced degrees of freedom mean the correlations are no longer significant at the 95% level for the annual mean cases at depths of 2022 m and 3800 m ( Fig. 4e and f). At the deepest depths (3800 m) we find poor correlations in both unfiltered and annual mean cases. Taken overall however, the correlations are encouraging, suggesting that the array estimates the maximum AMOC accurate to better than 0.5 Sv on annual timescales. Thus observed features such as the decline in the AMOC since ∼2005 and the major reductions in 2009 and 2010 are comfortably above the structural error (see Section 4.6 for further discussion of this point).
Beginning at the surface and working down the water column we see that in the top 100 m (Fig. 5a)   (see above). The bias changes sign at about 40 m depth, and occurs because the Ekman transport in the RAPID method is assumed to be depth independent in the top 100 m. The abrupt jump at ∼97 m again occurs as a result of the assumption of a constant Ekman-layer velocity. We note that in Fig. 5b the net ageostrophic transport bias integrates to close to zero at ∼100 m depth indicating that the net model ageostrophic transport in the Ekman layer is in close agreement with the Ekman transport calculated from the zonal wind stress (which is used for the proxy transport estimate). The other significant term in the top 100 m is the transport in the missing triangle at the eastern boundary which contributes ∼0.01 Sv/m on average (green line).
The unsampled regions remain unimportant between depths of 200-2500 m and the main contributing terms are the geostrophic reference level and compensation terms, roughly equal and opposite above 1000 m (Fig. 5c, red and cyan lines). Surprisingly, the ageostrophic term plays quite a significant role between 200 and 1000 m, but it is small everywhere below 1000 m. Between 3000 and 4320 m depth, the dominant term is the unsampled regions west of the MAR and at the eastern boundary. The compensation and reference level terms remain significant, but cancel each other. The reference level term results in a negative transport bias.
Below 4320 m the AABW, before mass compensation, is specified as the true mean transport over years 1990-2004 of the model simulation, and hence the mean bias over the 15 year period is composed solely of the mass compensation term. Fig. 5b and d show the terms integrated with depth, and indicate that the total bias in the Ekman layer peaks at about −1Sv (Fig. 5b), mostly caused by the assumption of a constant Ekman velocity. Over the whole water column, the maximum bias occurs at ∼3000 m depth and attains a magnitude of about −2.7 Sv (Fig. 5d). The AABW is overestimated by about 1.2 Sv at ∼5000 m depth. The reference level and compensation terms are of opposite sign and attain magnitudes of −12.2 and +12.2 Sv respectively at the deepest depths. Again, somewhat surprisingly, the ageostrophic flow is non negligible, peaking at −1.0 Sv at ∼700 m depth.
A useful way to assess the importance of each contributing term is to split the compensation term into portions which individually compensate the other terms plus a residual due to the fact that there is a net transport through the section (Section 3.3, Eq. (19)) as shown in Fig. 6. The effect in the upper 100 m is simply to cancel out the reference level term (Fig. 6a, b versus Fig. 5a, b). However at deeper depths (Fig. 6c, d versus Fig. 5c, d) the effect is to reduce the importance of the reference level term whilst increasing the importance of the unsampled regions term. As each term is individually compensated this gives an idea of the error associated with each term in the absence of the others. Thus if all the unsampled regions were instrumented then the overall streamfunction bias at 3000 m might be reduced by ∼1.1 Sv, but a reference level error of −2.2 Sv would remain. Note that the residual term is the compensation associated with the non-zero mean net transport through the section. For reference, summary statistics showing the importance of each term at four different depths are presented in Table 1.

Decomposition by region
The reference level bias term, T p -T G (red lines in Fig. 5a, c), can be decomposed into terms corresponding to the velocity along the lower boundary of the rectangles defined by the moorings as discussed in Section 3.4 (Fig. 7a, Eqs. (24) and (25). In this case it is easier to work upwards from the deepest depths. Below the depth of the top of the MAR (3000-4300 m), the bias is composed of two terms corresponding to the western and eastern basins respectively (solid red and blue lines). The deep southwards currents of the western basin combine with the rather weak currents of the eastern basin, to give a net negative transport bias. As we clear the top of the MAR two further terms contribute, the negative velocities (solid green) at the level of the top of the MAR and positive velocities (cyan) at the eastern boundary (the physical origin of the deep northwards flow at the eastern boundary is not clear). The MAR velocities are stronger than the northwards velocities at the eastern boundary and so the two new terms increase the overall bias at depths of 2000-3000 m. Of the remaining terms only that associated with the shallowest eastern boundary mooring contributes significantly between 2000 m and the surface (dashed green). The shallow southwards currents at this mooring make a surprisingly large contribution to the bias.
The direct impact of the unsampled regions (green line in Fig. 5a, c) is decomposed according to Eqs. (26)-(29) and the results are shown in Fig. 7b. The main extra information gained is in the lower levels where between 2500 and 4500 m the deep eastern boundary triangle contributes a strong positive transport, but the eastern MAR triangle makes a negligible contribution. The western MAR triangle exhibits a reversal in current direction at about 3400 m, reducing the positive bias above 3400 m but increasing it below this depth. Since at these depths there is bias compensation between the reference level term and the unsampled region term, correcting either one of these biases without correcting the other would actually increase the bias in the proxy streamfunction.
The other substantial contribution is from the ageostrophic term (Ivchenko et al., 2011). The term itself (blue line in Fig. 5a, c) includes the wind-derived Ekman transport minus the ageostrophic transport. Fig. 8a, b split the ageostrophic transport into contributions from the areas between each mooring pair (in this case it is more instructive to consider only the true model ageostrophic transport, T G , not the wind derived, depth constant, proxy Ekman transport, T EK ). In the Ekman layer (Fig. 8a) the net transport (solid black) is essentially the true, depth-varying, model Ekman transport, and the areas between the short mooring pairs at the eastern boundary (dashed lines) do not contribute much of this transport (solid black) across the section. The areas between the other four mooring pairs contribute to the Ekman transport in roughly equal measure, but interestingly the vertical profile of the Ekman transport changes character between the western and eastern basins: in the west the Ekman transport decreases monotonically with depth (solid red and green lines), whereas in the east it grows to a maximum value at around 10 m depth before gradually decreasing to zero (blue and cyan lines). A gradual deepening of the Ekman layer eastwards along the section is also observed in the model (not shown).
As already mentioned, the ageostrophic term below the Ekman layer (Fig. 8b), although small, is not negligible and contributes nearly 1 Sv to the overall bias. The largest values, peaking at 0.002 Sv/m occur in the top 500 m and much smaller values of <0.0004 Sv/m occur at deeper depths. The most interesting result is that the ageostrophic transport in the top 1000 m is spread roughly equally across the ocean regions, excepting the areas between the three shallow mooring pairs at the eastern boundary, and resides in the main thermocline.
The largest values of the ageostrophic velocity are found at the western boundary and in the Ekman Layer (Fig. 8c). At the western boundary, ageostrophic velocities of up to ±1 cm/s occur and much larger values are seen (as expected) in the Ekman layer. At the western boundary, the ageostrophic currents are a factor of 10% or less of the geostrophic currents, and generally of the opposite sign. To investigate whether the ageostrophic currents at the model western boundary are of a realistic order of magnitude, we take advantage of the fact that current meters measure the real-world transport between X 1 and a point known as WB3 at 76.5°W, 27 km east of the WBW used in the RAPID calculation. Geostrophic real world transport can also be calculated between X 1 and a full depth mooring at WB3. The observed geostrophic transports are calculated relative to the deepest common level and a reference level velocity is derived from the direct current meter measurements. Fig. 8d compares the directly measured real-world transport, including ageostrophic transports (thick blue), and geostrophically estimated real-world transports (thick black) with model equivalents (thin blue and black lines). The observed geostrophic transport and the full transport agree well. The geostrophic transport results in 0.4 Sv more northward transport shallower than 1100 m and 0.4 Sv more southward transport between 1100 m and 3400 m (3400 m is the shallowest depth used in the model). The model geostrophic transport results in 0.1 Sv more northward transport than the full velocities shallower than 1100 m and 0.33 Sv more southward transport between 1100 m and 3400 m. An over-strong ageostrophic flow near the western boundary was cited by Stepanov et al. (2016) as the primary reason for a discrepancy between model and RAPID estimates of meridional heat transport. However, the results quoted above indicate that the 1/12°resolution NEMO model used in this study does a reasonable job in capturing the magnitude and sign of the ageostrophic flow at the western boundary. Fig. 8d does, however suggest that the model deep western boundary current is too strong compared to the RAPID observations at these longitudes. For example the maximum southward transport in the model at this location is 0.0045 Sv m −1 at 2000 m depth. The estimate based on the observations suggests 0.003 Sv m −1 . Fig. 8e explicitly shows the model ageostrophic transport per unit depth as a function of depth and its standard deviation. The time variation of the ageostrophic transport at the western boundary is seen to be relatively small, in contrast with the results reported by Stepanov et al. (2016).

Time variation of biases
Finally, we examine the evolution of the bias in time. Fig. 9 shows the net bias as a function of depth and time. As expected from the plots of the mean bias, there is in general a dipole bias in the surface 200 m (negative at the surface and positive at deeper depths) associated with the representation of the Ekman transport as depth independent in the top 100 m (not shown, Fig. 9 shows the depth range 100-4389 m, as otherwise the Ekman layer would dominate in transport per unit depth, though not in net transport). Between 200 m and approximately 3000 m depth the bias is predominantly negative (i.e. the proxy underestimates the true transport) and below 3000 m it is predominantly positive with maximum amplitude at ∼300 m and 3700 m depth respectively (c.f. Fig. 6c). These compensating biases in depth may have implications for associated heat transport biases. The bias shows an interesting time variation. The bias is relatively low during the model spinup period, years ∼1978-1983 and subsequently displays marked interannual variations (2-3 year timescale) including large negative biases at ∼300 m in 1998, 2002 and 2004. Peak positive values occur at around 3700 m depth in 1998. After peaking in the period 1996-2005, the biases reduce considerably in the post 2006 period. We note here that Stepanov et al. (2016) showed that biases at depths greater than 500 m due to the RAPID methodology do not have implications for associated meridional heat transport biases, whilst the biases in the top ∼500 m are significant. Stepanov et al. (2016) reported that these biases were mainly due to the variability of the recirculation of the subtropical gyre (and eddy variability) at the western boundary. This may also be the case in our model, although we have not performed the analysis required to demonstrate this.  6. (a) Individually compensated biases in transport per unit depth (Sv/m) in the top 120 m. Net bias (black) and components associated with the reference level bias (red), unsampled regions (green), ageostrophic transport (blue) and residual (magenta) averaged over years [1990][1991][1992][1993][1994][1995][1996][1997][1998][1999][2000][2001][2002][2003][2004]. Note that the residual is not compensated. (b) streamfunction bias (Sv) in the same quantities (i.e. previous lines integrated w.r.t. depth), (c) as (a) for depths between 120 m and 6000 m, (d) as (b) for depths between 120 m and 6000 m.  Fig. 10 shows time series of the bias contributions from the reference level correction, ageostrophic transport, and unsampled regions, individually compensated as in Fig. 6, for different depths. In the upper layer (0-662 m) the variability is determined mainly by the reference level and unsampled region terms, with the reference level correction dominating on interannual timescales (Fig. 10a). There is a noticeable anticorrelation between the unsampled regions and the reference level term on timescales of ∼5 yr and higher; this is not unexpected as the two terms are both linked to the velocity in the missing regions. The ageostrophic term displays relatively small interannual variability. The same general picture is found at deeper depths (662-2865 m and 2865-4389 m, Fig. 10b, c), however the amplitude of both the net bias, and the interannual variability in the reference level term is much larger and there is a decadal timescale spinup effect particularly marked in the unsampled region bias (Fig. 10c).
Structural errors such as those investigated in this paper are very likely to have an impact on estimates of salt and heat transports in addition to mass transport. However the estimation of heat and salt transports is more complicated than the estimation of volume transports and a separate study would be required to quantify the biases in these quantities.
Nonetheless an order of magnitude estimate for the size of the overall bias in estimates of overturning heat transport ( MHT Δ ) can be provided by multiplying the proxy transport per unit depth, ∂ ∂ ψ z / p , by the zonal average of the temperature, θ, integrating with respect to depth, and comparing with a corresponding calculation using the true model streamfunction, ∂ ∂ ψ z / : where Cp is the specific heat capacity of sea water. The third term in the square brackets on the RHS of Eq. (30) is a constant velocity selected to compensate for the fact that the true streamfunction (unlike the proxy) yields a non-zero depth integrated transport. For our simulations, MHT Δ is negative and has a mean value of −0.16PW, a significant proportion of the estimated ∼1.2PW of heat thought to be transported northwards by the ocean circulation at 26°N (Johns et al., 2011), however this should be taken as indicative only: an extension of the present work is required to rigorously evaluate the structural error in RAPID estimates of MHT and its components due to missing regions, ageostrophy and reference level assumptions. Haines et al. (2013) identified that lack of measurements in the surface 50 m would effectively reduce the heat transport estimate by ∼0.1PW, which led to subsequent improvements in the surface extrapolation methods in RAPID to eliminate this bias.

Summary
We have performed a thorough analysis of the possible biases arising in estimates of basin-wide volume transport using mooring arrays such as RAPID. Three main sources of error are identified: those due to an assumed fixed reference level or level of no motion, those due to ageostrophic flows and those due to unsampled regions (e.g. Fig. 6c). Two other sources of error are found to make minor contributions. These are the assumption of a constant Ekman velocity (see Section 4.2) and the assumption of zero net flow through the section (see Section 4.1). The assumption of a constant Coriolis parameter between mooring pairs makes a negligible contribution.
The model predicts that the RAPID AMOC estimate (i.e. the maximum value of the streamfunction, which occurs at ∼900 m in the model) is likely to be too low by O(1-2 Sv). However our results also show that the standard deviation of the bias at this depth is small, O (0.3 Sv) on annual and longer timescales, compared to variability in the AMOC of O(2Sv), showing that the RAPID array is well suited to the studies so far conducted (e.g. Smeed et al., 2014McCarthy et al., 2012Fig. 7. (a) Decomposition of reference level bias in transport per unit depth (Sv/m) between mooring pairs averaged over years 1990-2004, western basin, X 2 -X 1 (solid red), eastern basin, X 4 -X 3 (blue), Mid Atlantic Ridge, X 3 -X 2 (solid green), eastern boundary, X 5 -X 4 (cyan), eastern boundary, X 6 -X 5 (black dashed), eastern boundary, X 7 -X 6 (red dashed), eastern boundary, X 8 -X 7 (green dashed). Locations  Duchez et al., 2016;Moat et al., 2016). We note that the bias can vary on interannual and decadal timescales, dependent largely on changes to the velocities along the baselines of the array mooring pairs and on the transports in the unsampled regions, hence the bias during individual events may be slightly larger or smaller than the mean bias. For example during the 2009-10 downturn, our simulation suggests that the bias was smaller than average (order −1.2 Sv at ∼660 m, Fig. 10a), and much smaller than the downturn itself (approximately 6 Sv lower than the long term mean for the period 1 April 2009-31 March 2010, McCarthy et al., 2012. However a scatter plot of bias versus AMOC at ∼900 m depth (the level of the maximum) shows no correlation between the AMOC strength and the bias, so there is no suggestion that the proxy bias is systematically larger during extreme events (Fig. 4).
The ageostrophic term in the model is comparable to corresponding estimates from observations in the western boundary (Fig. 8), however its size in the main thermocline is somewhat surprising. It would be of interest to further investigate this result by looking at the origin of this ageostrophy (e.g. nonlinear terms, viscosity) in the NEMO model and also quantifying the ageostrophic term in other models. The quantitative agreement at the western boundary between the model and observations is, however, encouraging (Fig. 8d). It should be noted that in the present study we were unable to evaluate the contribution to the ageostrophic velocity from correlations between changing thicknesses of the model grid boxes and modeled velocities and pressures. Although this is expected to be small, it is at present difficult to quantify precisely. We have excluded grid points adjacent to model boundaries in evaluating the ageostrophic velocity.

Historical context
Of the three main assumptions tested, each has been investigated before. For example Chidichimo et al. (2010) speculated on the difference between an eastern boundary tall mooring and instruments on the slope -equivalent to testing for a bottom triangle, and Ganachaud (2003) evaluated the uncertainty in basin-wide one-time hydrographic surveys due to unsampled regions, suggesting that this was negligible for those surveys. In contrast to Ganachaud (2003) our results find significant transport in the unsampled regions at depth at the eastern boundary and also west of the mid-Atlantic Ridge. To our knowledge, this has not been done systematically for the remaining uninstrumented regions of the current RAPID array. We were not able to test for the influence of the unsampled region at depth in the western boundary  [1990][1991][1992][1993][1994][1995][1996][1997][1998][1999][2000][2001][2002][2003][2004]. Full section, X 1 -X 8 (black), western basin, X 2 -X 1 (red), eastern basin, X 4 -X 3 (blue), Mid Atlantic Ridge, X 3 -X 2 (green), eastern boundary, X 5 -X 4 (cyan), eastern boundary, X 6 -X 5 (black dashed), eastern boundary, X 7 -X 6 (red dashed), eastern boundary, X 8 -X 7 (green dashed). Locations, X 1 , X 2 …X 8 along the RAPID section are indicated on Figs. 2a and 2b), (b) as (a) for depths between 120 m and 6000 m, (c) modeled ageostrophic velocity (shaded, m s −1 ) along the RAPID section close to the western boundary. Contours of geostrophic velocity (m s −1 ) are overlaid, dashed contours representing southwards velocities. The thick vertical dashed line represents the eastern limit of the western boundary wedge. (d) Observed net transport per unit depth in Sv/m (thick blue) and observed geostrophic component (thick black) and model equivalents (thin blue and black respectively). (e) Model ageostrophic transport per unit depth (solid red) and its standard deviation (dashed red) both in Sv/m. Observed transports are evaluated over the longitude range 76.74-76.5°W averaged over the period 2004-2014, modeled transports are evaluated over the longitude range 76.67-76.42°W averaged over years 1990-2004. due to the boundary between NADW and AABW being located at too shallow a depth in our model simulation.
The reference level assumption has been discussed many times and at length. In particular the work of Baehr et al. (2009)  B. Sinha et al. Progress in Oceanography 160 (2018)   than that of Baehr et al. (2009) and brings out the underlying reason for the sensitivity to the reference level, the precise calculation of the bias in terms of currents at the bases of the mooring pairs, and hence the relative importance (and potential error compensation) of particular mooring pairs, and the use of the absolute geostrophic velocity determined using the model sea surface height field. We emphasize that no other study has capitalized on the availability of the sea surface height in a model to construct the absolute pressure gradient for comparison with geostrophic calculations (the theory has of course been developed and used in other contexts, e.g. Hughes et al. (2013) and Willis (2010) for the 40°N AMOC estimates, and Ganachaud (2003) used the SSH to provide estimates of the variability of the velocities in the real-world reference level velocities). Our analysis relates the bias to the velocity between particular mooring pairs. This diagnostic, i.e. the velocity between the moorings (e.g. V21 and V43 in Fig. 2a) varies with model drift and shows how the reference level bias changes with time (e.g. Fig. 10a). Hirschi et al. (2003) tested the validity of the RAPID method in the presence of a drift, but their purpose was different. There it was to identify if trends can be identified by the RAPID methodology (they can), here it is to quantify how the bias changes as a result either of real trends, or as an artifact of model spinup.
Although many of our results may not seem very surprising, this study needed to be done in order to estimate the importance of these biases in a direct analogue of the actual RAPID array. High resolution studies performed prior to the deployment of the real RAPID array were intended as a proofs of concept (Hirschi et al. 2003, Baehr et al. 2004 and used idealized array configurations. The real array was subsequently deployed in a configuration different to these earlier studies. More recent high resolution studies (e.g. Roberts et al. 2013;Duchez et al. 2014) have focused on studying how different methodologies to measure the AMOC are likely to perform rather than on the structural error inherent to the actual mooring setup in RAPID. Hence, despite the large body of literature dealing with monitoring methodologies the question remained: Can we quantify the biases in the real array?

Importance of horizontal resolution
Our higher resolution model suggests a larger role for ageostrophic transport in the main thermocline than expected from the Baehr et al. (2009) analysis which used 1°resolution models. Lower resolution models display very weak nonlinearity and in addition lower resolution in the horizontal may result in weaker vertical gradients and lower viscous drag. Both processes can lead to ageostrophic transports and may therefore be stronger in the higher resolution model. An additional effect of higher resolution is the ability to resolve stronger, narrow boundary currents (including along topography such as the MAR). These fall either in the unsampled regions (Fig. 7b), or at the bases of the mooring pairs, contributing to the reference level bias (Fig. 7a) and hence there is potentially a modified balance between (and a modified net effect of) these two terms compared to lower resolution models. A further effect of higher resolution is that the Antilles Current is more closely confined to the boundary and therefore better captured by the WBW transport. In lower resolution models the Antilles Current is more diffuse, hence a larger fraction of its transport is captured by the part of the array east of the WBW, which relies on geostrophy and is subject to the reference level error. On all these counts we expect different contributions to the bias in a high resolution model compared to lower resolution.

Implications for ocean modeling
A key point is how far we can trust the model estimates of the bias terms. The major bias terms originate essentially as circulation features in the unsampled regions of the model. These regions are all along the boundaries of the section: east, west and bottom. It is as yet unclear how well the model captures the real circulation and its variability in these regions. Hence both observational and intermodel studies of the circulation in these regions would be effective in reducing uncertainties.
The spinup period of the model introduces spurious drifts and associated biases, and the analysis highlights potential shortcomings of the model. Chief among these is that since the proxy transport underestimates the actual transport in the model, then the mean transport in the real North Atlantic at 26.5°N is likely 1-2 Sv above the estimate from the real RAPID array (17 Sv). Hence the model estimate of the AMOC at 26.5°N (∼15.2 Sv, Fig. 3) may underestimate the real transport by 3-4 Sv.
To our knowledge, a detailed and rigorous study such as this one has not previously been performed and given the importance of the RAPID array to studies of Atlantic and global climate variability, there is a need to gain in-depth understanding of possible error sources in the AMOC observations arising from the current observational setup. Since such work would be too onerous in the real ocean, state of the art ocean models provide an ideal testbedjust as they did in the early phase of the RAPID project when the original observational strategy was being devised. Earlier studies have shown that the model used here realistically simulates the key features of the North Atlantic ocean circulation and therefore we are confident that our findings largely translate to the real ocean. Nevertheless it would be highly desirableand instructiveto perform the analysis on a variety of high-resolution (eddy resolving) models as factors such as grid-type, coordinate type, resolution, surface forcing and parameter choices are likely to subtly affect the relative contribution of different processes to structural error estimates of the AMOC. We note that Stepanov et al. (2016) have performed a similar study using a higher resolution model. Whilst their paper is more focused on the differences between the model and the real array, their equivalent of our Fig. 3a indicates that in their model, the RAPID proxy appears to be ∼3 Sv higher than the true AMOC whereas in our model the RAPID proxy is less than the true AMOC. It is unclear why this is so, but our framework applied to their model would clearly bring out the reasons.

Potential improvements to the RAPID array
Better characterization of the mooring baseline velocities and sampling of the regions currently uninstrumented would improve the RAPID transport estimates, although error compensation needs to be guarded against. Large reduction in one of these biases without corresponding reduction in the other could lead to reduced accuracy. Some historical epochs appear to have much larger biases than others, suggesting that the accuracy of the RAPID method might be affected by some physical processes that vary in time (e.g. Stepanov et al., 2016): these could potentially be corrected for, or a time varying accuracy could be quoted. Use of correlations between the various sources of error and the measurements themselves could also be utilized to reduce error (for example using correlations between the geostrophic transport and the unsampled regions to extrapolate into the uninstrumented areas).
Given the importance of the RAPID array and the development of basin-wide monitoring arrays at other locations (Lozier et al., 2017;Ansorge et al., 2014) the community needs a framework for modeling structural biases in these arrays based measurementsthis is the underlying philosophy of this paper. In this way advanced ocean models can be integrated into the design and redesign of observational arrays. And to emphasize it one more time, the exploitation of the absolute geostrophic velocity to bring out the reference level bias has not, to our knowledge been applied to basin-wide array measurement before.

Concluding remarks
In a generic way, error analysis is a characteristic and important B. Sinha et al. Progress in Oceanography 160 (2018) 101-123 part of science, if perhaps less glamorous than other activities. More specifically, accurate estimation of the biases are important as the array transports are being used to calculate transports of heat, freshwater and carbon. If basin-wide budgets of these quantities are to be computed accurately, then the mean bias and variability in these quantities need to be well constrained.
The sources of error considered in this paper do not compromise the principle of the RAPID array and other arrays like it. Indeed more accurate characterization of the errors associated with basin-wide measurements of mass, heat and salt transports will enhance their usefulness and impact on scientific research and wider society.
Finally, our purpose in this paper is not to question the methodology or robustness of the RAPID array (that has been demonstrated on many occasions previously, and the RAPID array has emerged unscathed from a number of apparently serious criticisms). Instead, our purpose is to rigorously quantify (using a well-respected ocean model and a perfect model approach) the structural biases inherent in the method and their sources and relative magnitudes. This in no way compromises any of the conclusions of previous studies, and in fact emphasizes their validity.
Although many studies have been performed on single sources of bias, this is the first attempt to draw them all together into a unified framework using an eddy resolving ocean model. Our work is unlikely to be the last word on the subjecteven higher resolution models are starting to be used for this purpose (e.g. Stepanov et al., 2016), but in order for these to be useful, we need to have an agreed procedure for systematic comparison of models-model and model-observation differences. and = allows us to write The first term on the RHS represents the geostrophic transport between X w and X e , whilst the second represents the ageostrophic transport due to the variation of the Coriolis parameter. Note that the geostrophic term is independent of the absolute value of pressure (P e and P w ) and uses an average value of the reciprocal of the Coriolis parameter. It is thus important not to use the formula − ( ) Appendix B. Bias due to assumption of a fixed reference level for dynamic height calculations We first consider an idealized case with single moorings on the eastern and western boundaries and one on each side of the MAR (see Fig. B.1). Dynamic height, D, relative to a reference level Z REF is given by where g is the acceleration due to gravity, ρ is the in situ density and ρ 0 is a reference density. Then the pressure, P, at a depth, z, is given by Now, if we consider the absolute geostrophic transport, T, in the west basin area bounded by = x X X , a b and = z z 0, REF , then where the absolute geostrophic transport per unit depth, (multiplied by the Coriolis parameter, f), is given in terms of the SSH, h, difference plus dynamic height difference (subscripts refer to the mooring positions defined in Fig. B.1), and hence as the absolute geostrophic velocity (v ab ) at the base of the mooring pair times the distance between the moorings.
Above the level of the MAR, the proxy geostrophic transport, T proxy is derived relative to a fixed reference level using Eq. (7): The third and fourth lines represent the absolute geostrophic transport: Below the level of the MAR and similar analysis yields For our second case we consider a situation with a single mooring to the west and a stepped profile (N steps) at the eastern boundary with no MAR (see Figs. B.2 and B.3).
The proxy transport above the N th step is given by The terms in h can be taken outside the summation to give  Combining the two cases, we see that the reference level error at any depth is composed of the sum of the transports from the geostrophic velocities at the baselines of the mooring pairs below that level. This is illustrated graphically in Fig. B.1.
The preceding analysis brings out the difference between the true transport and one calculated using geostrophy between the mooring pairs, however this is not the error in the RAPID array estimate because the compensation velocity, c, introduced in Section 3.1 partially offsets it. Thus the actual error between moorings n and n + 1 is proportional to (v nn+1 -c). Since c is constant but v nn+1 varies with n, this underlines the fact that the compensation cannot perfectly cancel out the reference level error.  Illustration of an example domain with tall eastern boundary mooring at X 0, spanning the entire water column from H 0 to H N and short moorings at the stepped eastern boundary at X 1 , X 2 …, X N . The mooring at X 1 spans the water column from H 0 to H 1 and the sea surface height at this location is h 1 and so on for the moorings at X 2 , X 3 …, X N . below the MAR (levels 67-59 inclusive), where ∼ f is the average value of the Coriolis parameter and ρ 0 is a reference density (values for these are given in Appendix C.4).
The model equation of state is based on the UNESCO equation of state (UNESCO, 1983further details of its implementation are given by Madec (2008)) and calculated as a function of model prognostic variables salinity and potential temperature and depth (assumed fixed for the virtual moorings, but in fact the z levels do varysee Appendix C.4) C.3. Calculation of geostrophic velocity based on level of no motion In order to calculate a RAPID-style proxy meridional streamfunction from the model it is necessary to specify the timeseries of Florida current transport, the flow at the western boundary wedge and the AABW. These are calculated as follows: Firstly the meridional velocities corresponding to model gridpoints i = 2471-3301 and j = 1823-1824 and k = 1-75 were extracted for all 5-day means between January 1978 and December 2010. The velocities were then averaged onto the locations of the zonal velocities corresponding to i = 2471-3300 and j = 1824 (this is explained further in Appendix C.5).
For the combined Florida Strait and WBW transports, these averaged velocities were then summed (area weighted) over the range i = 2471-2525, j = 1824, k = 1-75. For the AABW the velocities were summed (area weighted) over the range i = 2471-2525, j = 1824, k = 68-75 and averaged over the years 1990-2004 to avoid the model spinup period. The area weighting used for the sums varied with time because of the time dependence of the gridbox vertical thicknesses.
A constant value of the Coriolis parameter, ∼ f , equal to 6.47 × 10 −5 s −1 was used, obtained by averaging the Coriolis parameter at each gridpoint for the range i = 2471-3300, j = 2823-1824. The reference level density, ρ 0 is taken as 1035 kg m −3 and the acceleration due to gravity, g, was given a value of 9.81 m s −2 .
The Ekman transport at each gridpoint was based on the 5-day mean model zonal wind stresses extracted from the model output for the range i = 2526-3300, j = 1824, i.e. excluding the Florida Strait and the Western Boundary Wedge.  The moorings are located so as to avoid partial gridpoints in the calculation of the proxy transport (Appendix C.1) and absolute pressure gradient and the ageostrophic velocity (Appendix C.6), however they must be taken into account in the calculation of the Florida Strait, WBW, AABW and net transports and in the calculation of the transports in the unsampled regions (Appendix C.3).

C.5. Calculation of reference level bias terms
This is the discrete analogue to the result from Appendix B. Thus for the reference level bias associated with moorings 7 and 8, the contribution is This is the pressure difference (divided by g) times the distance between the moorings (X 8 -X 7 ) and must be multiplied by the depth (and by 1/ ∼ ρ f 0 ) under consideration to give a transport bias, e.g. at a depth of ∼700 m (k = 42) where the bottom of mooring 8 is located, we must multiply by ∼700 m. Similar calculations must be performed for each mooring pair and the contributions summed up. Note that at the surface all mooring pairs will make a contribution, but at deeper depths, only those mooring pairs which reach below the depth under consideration will contribute. At the reference level, Z ref , only mooring pairs X 1 , X 2 and X 3 , X 4 will contribute.

C.6. Calculation of ageostrophic velocity
The ageostrophic velocity is calculated as the difference between the model prognostic velocity and the model pressure gradient. As it is a small difference between two large numbers, it is essential that when calculated offline as in our study, the prognostic velocity and pressure gradient are calculated in exactly the same way as the model solves the equations of motion.
NEMO utilizes a C-grid which means the positions of the discretized variables relative to each other are as in Fig. C.6 (after Madec, 2008). In the prognostic equation for the meridional velocity at gridpoint (i, j), u i,j , the pressure gradient is obtained by taking the difference between the pressure at the adjacent gridpoints, P i,+1j -P i,j . The meridional velocity at the same point is given by 0.25 * (v i,j + v i,j-1 + v i,+1j + v i+1,j-1 ) and the Coriolis parameter as 0.25 * (f i,j + f i,j-1 ). In fact the model itself uses a slightly more complicated formula for the product fv at gridpoint (i, j), involving the average of four terms such as v i,j (f i,-j + f i,j + f i,j-1 )/3. This is necessary in the model to ensure exact conservation of vorticity and enstrophy. At 26.5°N, the model gridlines follow latitude and latitude lines very closely and there is little variation in Coriolis parameter from gridpoint to gridpoint, hence our slight simplification introduces negligible error to the estimation of the ageostrophic velocity.

C.7. Model parameter settings
These are listed in Table C.7.1.