ACCURATE SENTINEL-2 INTER-BAND TIME DELAYS

: Accurate knowledge on the time-delays between detector bands is paramount for the observation of dynamic features in Sentinel-2 imagery. Sentinel-2 inter-band time delays have been characterised analytically and empirically on Level-1B images with help of the rigorous geometric model. Two major contributors have been identified: optical distortion and satellite altitude inducing across- track and along-track delay variations up to +/-2% with respect to ESA handbook fixed value. A harmonic attitude perturbation induces a 0.15% peak to peak temporal modulation on delays. An original method is proposed in order to estimate inter-band delays from L1C orthorectified and tiled products. The method relies on the orbital ephemeris embedded in L1C metadata and a lookup table. The performance assessment exhibits a maximum error of 1% and a rms error of 0.28%.


INTRODUCTION
In recent years, large scale, systematic, Earth Observation programmes like the Sentinel series or alike, have become the new norm to observe and monitor land and ocean processes, climate impacts, and anthropogenic changes (European Commission, 2015). The systematic nature of these image acquisitions that last for several years to decades, enables investigation of time-series with unprecedented temporal resolution. While the Sentinel-2 satellite programme was initially designed to monitor time-varying processes instantaneously at a given moment (meaning a static observation that might change in time), it was not foreseen that the Sentinel-2 satellites could be used to observe instantaneous but dynamic physics such as propagating features (Kudryavtsev et al., 2017;Yurovskaya et al., 2019;Bergsma et al., 2019).
Exploiting a small, but sufficient, inherent time-lag in pushbroom sensors between image bands, enables the detection of propagating features, both natural like ocean currents (Yurovskaya et al., 2019), waves (Kudryavtsev et al., 2017), wave-derived bathymetries (Bergsma et al., , 2021, clouds, volcano plumes (de Michele et al., 2019) and unnatural objects like inflight airplanes (Liu et al., 2020;Heiselberg, 2019). The dynamic property one is after in these cases is a velocity or celerity in some sort. Considering that the time-lag is relatively small in comparison to the minimum pixel size and the target physics (in the case of ocean currents, wave and cloud propagation), any error made in space and time is a first order linear error, and requires hence to be measured accurately, (and precisely). The spatial part is often methodology, and/or process depends, hence, in this work we focus on the common factor: the time-lag between image bands.
The time-lag between image bands is stated (in absolute form) in the ESA manual for Sentinel-2 (MSI instrument manual) but considering the complex focal plane configuration of the Sentinel-2 instrument, consisting of multiple detectors, these given time-lags should be considered with caution. If not processed properly, the given, constant, time-lag between image bands results in different derived velocity components per detector band, that are not necessarily just a sign difference. A time-lag between the image bands per detector is therefore required to overcome this issue. At the same time, time information at this level of detail is not delivered, nor can be derived directly from the standard available Sentinel-2 products, distributed by ESA: L1C, L2 and higher. Consequently, in most of former cited works, inter-band time lags are taken from ESA table. A very few works worked about a finer estimation of Sentinel-2 time delays. (Yurovskaya et al., 2019) estimate Sentinel-2 delays thanks to product metadata such as satellite altitude and speed ephemeris, as well as azimuth and zenithal viewing angles for each spectral band. The authors claim a 1-2% accuracy, but did not detailed the assessment, nor the methodology. (Liu et al., 2020) proposed an automated method to detect flying aircraft through the use of the inter-band time lags of Sentinel-2 measurements. This work relies on ESA table for inter-band time offsets, but they measured the relative time offset between adjacent detectors of B2 band thanks to images for which an airplane is sensed on two different detectors of the focal plane array, without knowledge of the aircraft speed and height.
In this work we propose a methodology, derived from L1B level products, to precisely compute a time-lag in between image bands and per detector for L1C and higher level products, with the meta-data available in the distribution of these products. The outline of the paper is as follows; the subsequent section describes the Sentinel mission, configuration and products, succeeded by a theoretical analysis of the time-lag calculation, including all constituents and potential error sources. We then demonstrate the empirical derivation of the time-lag using L1B level Sentinel-2 products, its dynamics over an orbit and the variability between detectors and the differences found depending on orbital position, altitude, and latitude. Using this knowledge, section 5, draws an approach to calculate the timelag per detector per image band pair for the distributed L1C or higher level products.

Sentinel-2 overview
The Sentinel-2 mission (Drusch, 2012) main-objective is to provide high-resolution optical imagery for the operational monitoring of land and coastal areas. Two identical satellites named hereafter S2A and S2B are maintained in the same sun synchronous 786 km mean orbit with a phase delay of 180° providing a revisit time of five days at the equator. Sentinel-2 satellites acquire images in a classical push-broom fashion with 12 linear detectors. The on-board telescope (MSI instrument) is designed with a staggered arrays focal plane. In order to acquire different spectral bands, the linear arrays are shifted in the telescope focal plane along the satellite along-track (ALT) direction (see Figure 1). Also, due to limitation of the size of the detectors, the 290 km field of view is obtained by mosaicking the detectors in the across-track (ACT) direction. Two focal planes optically conjugated co-exist: one for VNIR spectral bands, and one for SWIR spectral bands for a total of 13 bands. The VNIR spectral bands arrays (resp. SWIR arrays) lie on the same detector (chipset). Therefore, it is likely that the spectral band arrays are collinear in each detector. The spectral filters sort order is inversed between even and odd detectors. Though all detectors share the same telescope, the resulting images do not share the same ground resolution. Depending on the spectral band, 10m, 20m, or 60m images are distributed to users thanks to on-ground processing chain (Baillarin, 2012). Due to Time Delay Integration (TDI) technology, a yaw steering of the platform must be done so that ALT direction is perpendicular to the detector's arrays (see Figure 1). We call this orthogonality the "yaw steering condition". According to the Operational Concept Document (OCD) the main line of sight (LOS) of MSI is pointing towards the Earth' center. According to OCD,

Sentinel-2 products
Sentinel-2 images are de-spatialized into several levels of products. L1A and L1B products are dedicated to expertise while L1C and L2 products are delivered to end-users.
L1B images are radiometrically corrected and stored in native raw geometry, per detector. Therefore, 13x12 images are produced per L1B product on a datastrip. The L1B product metadata embeds the rigorous refined geometric model. Geometric refinement is performed through ground control points (Baillarin, 2012). Even without refinement, considering the very good absolute geolocation accuracy (of ~1 pixel), the absolute date of any image line of any spectral band can be retrieved with a very high accuracy (at least better than 2 ms depending on geometric refinement) and the relative dating accuracy is better than 0.15 ms (Languille, 2015).
L1C and L2 products provide georeferenced ortho-rectified multiband images with a sub-pixel interband registration accuracy for static object (Languille, 2015). L1C and L2 products are tiled in 110 km by 110 km tiles. A tile typically encompasses several detectors. A detector mask (MSK_DETFOO) allows to superimpose the detectors' footprint on the images. No detailed dating information specific to the tile and/or detector is available in the metadata.

THEORETICAL ANALYSIS
The following section deals with the analytical expression of the Sentinel-2 inter-band time-delays. The aim is to identify the main contributors to the time-delay and assess the expected variability over time or in the field of view.

Inter-band sensing
Let us consider a ground point M which has been sensed by spectral bands Bi at time tMi, i being the Sentinel-2 band number in intervals {[1-8], 8a, [9-12]}. Due to the focal plane layout and the yaw steering guidance law, this point has been sensed by the same detector number for all bands (see Figure 1). Two chipsets may be involved: one for the VNIR focal plane and one for the SWIR focal plane. For the sake of clarity and without loss of genericity, we will consider only one physical detector. We will also consider that the number of pixels is the same on all spectral bands. On this arbitrarily considered detector, the staggered spectral bands are spatially separated in the satellite track direction. At the exit pupil of the telescope there is an angular separation between any spectral bands couple which depends on the (Bi, Bj) couple involved. Consequently, the instantaneous ground footprints of the spectral bands arrays are also spatially separated in the satellite track direction. Taking into account the satellite ground track speed, point M has been sensed at different dates tMi and tMj by arrays Bi and Bj ( Figure 2).

Analytic formulation
Given the maximum time delay of 2,6s, the instantaneous spatial separation between spectral bands can be considered constant during this time interval. For the same reason the satellite velocity at times tMi and tMj can be considered equal. The platform attitude is geocentric-guided with null pitch and roll angles, and we suppose no attitude perturbation. The pinhole model is used to define optical projections. Let PM be the conjugate plane of the focal plane (Pfp) in which M lies illustrated in Figure 3. Let ⃗ be the ground velocity of the detector's projection at point M in an earth-fixed frame. Let ⃗ be the oriented orthogonal distance between the projected linear arrays Bi and Bj in PM. The time delay can be estimated from the following equation: The ground spacing ⃗ depends on the angular separation between Bi and Bj arrays in the focal plane, the telescope optical distortion, the satellite altitude, and the differential atmospheric refraction. Since atmospheric refraction is negligible with respect to incidence angles (Noerdlinger, 1999), the atmospheric refraction is not considered in the analysis here. ⃗ depends on the satellite apparent celerity ⃗ , and the satellite altitude Hsat* with respect to point M. The angle between ⃗ and ⃗ must be taken into account if the yaw steering condition is not met on point M. It is to be noted that yaw steering guidance cannot satisfy the yaw steering condition on every point of the swath if the telescope experiences optical distortion or if the detectors are not properly aligned.

Figure 4
Satellite apparent orbit and ground track. Projection of satellite velocity depends mainly on satellite altitude Hsat.
The quasi-circular orbit and the geocentric pointing constraint imposes a rotation of the satellite body such that the MSI main LOS follows the Earth' center ( Figure 4). Therefore, the rotation velocity ω of nadir ground track equals: , 2 where Vsat is the apparent satellite velocity, Rt is the local ellipsoid radius and Hsat is the satellite altitude w.r.t ellipsoid. The nadir ground track velocity at ellipsoid surface is then: where C is the Earth center and N the nadir point ( Figure 4). Considering now the velocity at point M, from Figure 3 one can deduct that neglecting ALT distance M from NP yields a good approximation of : Equation 5 shows that is always smaller or equal to nadir velocity. Indeed, Figure 4 illustrates that the edge of swath ground track perimeter is smaller than at nadir but the orbital period remains equal to the orbital period at nadir. By estimating the ground velocity ( ) for the Sentinel-2 mission, we find that the expected velocity range within the swath is -0.02% due to earth curvature. Note that taking mean earth radius Rm instead of Rt yields an error of same order of magnitude. From here onwards we will use the simplified equation to approximate the ground velocity: This allows us to express the orbital velocity on the ground but (6) does not take the yaw-steering into account while a yaw steering guidance imposes a yaw angular speed. The supplementary velocity imposed to M by yaw-steering is equal to: -.
where / is the yaw speed of the satellite illustrated in Figure 4. Figure 5 shows the yaw angle speed computed by an orbit simulator based on a Keplerian propagator. Figure 5 shows that the yaw speed is a quasi-linear function of the latitude in the range -80°/+80°. Its impact on ground speed is not negligible as will be illustrated on real data. Hence the total ground velocity ( ) is then a combination of the orbital velocity and yaw steering speed following:

Optical distortion
The main contributor to inter-band ground distance 23 is the angular separation between Bi and Bj sensing arrays in the object space (telescope exit pupil). This separation depends on the focal plane arrays assembly as well as on the optical distortion. It can be formulated as the line of sight (LOS) angular distance between Bi and Bj at the telescope exit pupil and it can be evaluated thanks to the in-orbit geometric calibration whose accuracy is better than 0.1 pixel (Dechoz, 2014, Languille, 2015. Indeed, the LOS of every pixel for every band and every detector is attached to L1B expertise products. Figure 7 shows the S2B LOS extracted on spectral bands B2 and B4 for every detectors and projected on across-track (Ψx) and along-track (Ψy) components. Nominally these quantities are defined in the MSI instrument reference frame but thanks to yaw steering guidance MSI axes are aligned with along-track (ALT) and across-track (ACT) axes. By definition, ALT direction is collinear to nadir ⃗ . Therefore, Ψy is the component which directly impacts the time delays.
Note that 45 6 7 , 45 6 / is proportional to the array's footprints projected on plane PM (Figure 3). Note in Figure 1 and Figure 7 that the order of staggered spectral bands depends on the detector number parity. Following the detector order, we expect a periodic sign change in the delays. This is why all delays will further be given as absolute values. The S2B LOS exhibits an optical distortion pattern which depends on the ACT position in the image. As optical design is identical between S2A and S2B, the same pattern is expected for S2A. The other consequence of optical distortion on time delays is that ⃗ and ⃗ vector cannot be perfectly aligned in the field of view. Indeed, LOS of detector #12 reaches 2.8° of tilt with respect to a distortion-free telescope. In order to evaluate the impacts on time delays, the ALT LOS difference is evaluated with the following formula for same pixel indexes in Bi and Bj LOS arrays: ∆ 456 / 89: 45 6 / ;< 89: ) 45 6 / ;= 89: , 9 where k refers to the detector array index. Referring to ⃗ definition, note that ∆ 456 / is not strictly equal to ALT component of ⃗ because 6 7_@2 89: and 6 7_@3 89: are not equal due to distortion in ACT direction, but since ∆ 456 7 is 10 times less than ∆ 456 / , the approximation is very good. Figure 8 exhibits the absolute ALT LOS difference A∆ 456 / A computed on spectral bands B2 and B4 for every detector. We conclude that for this contributor, we expect larger time delays at edge of the field than the middle of the field by a factor +3.25%.

Figure 8:
Absolute difference of B2 and B4 ALT LOS component. The swath range ratio due to distortion is 3.25%.

Satellite altitude
Referring to Figure 3, the inter-band array ground distance ⃗ projected ALT is deducted from the LOS angles: And since V HIJKLM ⃗ is ALT oriented, we can rewrite the delay equation (1): Since Hsat* and depend on satellite altitude Hsat, temporal delays have a strong dependency on the satellite' altitude. The satellite altitude (Hsat) modulation is due to: • Orbit eccentricity, • Ellipsoid flattening, • Surface relief.
Following Table 1, Sentinel-2 orbits have an eccentricity of 1.2E-3 and 7167 km for semi-major axis, which yields a quasi-circular orbit translated by 8km towards south pole with respect to Earth's center.
The WGS-84 flattening induces 21.4 km of altitude difference between equator and poles. These two contributors are the major ones. The combination of the two yields a 30 km range with a minimum altitude of 788 km for ~15° latitude and a maximum altitude of 818 km for the maximum latitude of 81°. The temporal variation is of very low frequency. Furthermore, Hsat* also varies across-track, as can be seen on Figure 6 as a static variation due to earth curvature.

Estimation of static ∆N variations
Among all the identified contributions, the static ACT contributors are: • Optical distortion • ACT altitude variation due to earth curvature From section 3.4, optical distortion contributes +3.25% to ∆ . Considering Sentinel-2 swath, from equation 11, earth curvature Hsat* contribution is +0.21%. From equation 6, we conclude that contribution is -0.02%. As can be seen in the simulation result presented in Figure 11, the total expected ACT range is +3.2%.

Estimation of dynamic ∆N variations
From equation 11, we find that the satellite altitude and ground track velocity are two dynamic contributions. We have listed in section 3.3 the Hsat contributors, being mostly the orbit and earth flattening. This 30km flattening yields an orbital range of 3.6% peak to peak but most users work on restricted +/-60° latitude, which reduces the altitude range to ~2%. Terrain height should be taken into account depending on the application. For bathymetry, there is obviously no relief contribution. From equation 8, orbital variations are due to and Hsat modulations. depends on earth rotation and satellite altitude. Indeed, apparent velocity inherits earth rotation being maximum near the equator. modulation due to Hsat follows Hsat profile, being maximum around +15° latitude. Simulations of nadir profile with Keplerian propagator yields an orbital range of 0.7% being maximum near the equator. The total expected ∆ range on Sentinel-2 orbit is 4.3% (3.5% between -60°/+60° latitudes).

Pointing errors
Sentinel-2 pointing error is specified to be maximum 2km from nadir (SMRD, 2012). Such pointing errors yield negligible delay variations.

Synthesis
The theoretical Sentinel-2 inter-band delays have been evaluated thanks to the knowledge of MSI LOS and the orbital parameters. All contributors have been identified and evaluated individually. Referring to the ESA time table which gives one delay per band couple, we have shown that inter-band ∆ undergoes a static ACT range of 3.2% as well as low frequency dynamic variations up to 4.2% (without relief contribution). The two main contributors are the instrument distortion and the ground-to-satellite altitude Hsat which impacts both the inter-array ground distance and the ground track velocity .

EMPIRICAL DELAYS MEASUREMENTS
The aim of this section is to check empirically on real Sentinel-2 images the theoretical statements of section 3 and to establish the true delays experienced by Sentinel-2 bands.

Principle
Keeping in mind that L1B expertise products are in native raw geometry and embed the rigorous physical geometric model, it is straightforward to compute ∆ on chosen ground targets.
For a given detector, given a pixel of index (ci,li) in image Bi, we use the geometric model to compute its ground coordinates on ellipsoid thanks to collinearity equation. Then the geometric model of image Bj is used in order to compute the coordinates (cj,lj) ( Figure 9). Since MSI is a pushbroom instrument, the acquisition dates of the two pixels in Bi and Bj are linearly related to their line index. Finally, the delay is simply the difference between these two dates. This process is applied 15 times for each of the 12 detectors: 5 line indexes distributed over the image length, and 3 columns distributed over the detector swath: first, middle, and last pixel. Delays have been computed for all band couples. Figure 9: Principle of delay estimation through image colocation illustrated on B2 and B4 bands.

L1B image set
The L1B image set is composed of 13 S2B products, and 2 S2A products. The ALT length of these products is variable, from 300km to 6500km. The acquisition dates range from 2017 to 2021 for S2B. S2A images are acquired in 2016 and 2017. LOS calibration files have been harmonized with the LOS version of 30 th october 2014 which is applicable to all products. The spatial repartition of the products is meant to be spread in latitude in order to evaluate the dependence on satellite altitude and latitude. The 840 measured delays are located in Figure 10.   Figure 11 shows the static ACT delays of band couple B2-B4 found for a particular date in a particular S2B product located in Spain. The same distortion pattern as Figure 7 is present. The total measured range is 3.2%. Figure 11 also presents results of a simulation of ACT variations, given empirical Hsat and values computed via the L1B product ephemeris. The good agreement validates our ACT static delay model. Simulating the true orbit would require an accurate orbit propagator which is out of the scope of the study.

Dynamic along-track empirical variations
The delays have been computed for the central pixel of detector #6 over the 13 S2B products for band couple B2-B4. The latitude range is -20°/+80°, which allows us to measure the range of delay due to Hsat and variations (blue dots in Figure 12). On Figure 12 is also plotted (orange dots) the Hsat / ratio which should be strictly proportional to the delay referring to equation 11. Hsat and have been estimated from the product ephemeris. Correlation is indeed perfect. The delay range is 1.9%, in line with our estimation on this range of latitude. It is to be noted that majority of the delay range should theoretically occur at southern hemisphere due to the orbit, but since the minimum latitude is -21° in our experiment, the range is much lower than the total expected 4%.   Figure 12 has been done with 60 points distributed over a high range of latitude. Figure 13 illustrates the same experiment with a high density of points on a shorter latitude range, and on only one L1B product. A periodic sine modulation is clearly visible, with a period of 30s and an amplitude of ~2ms. This frequency is the one observed during S2B in-orbit commissioning and identified as a solar panel vibration inducing a sine attitude perturbation at 0.032Hz (Languille, 2017). Indeed, a sine attitude pitch perturbation induce a sine ground speed modulation, hence a delay modulation. The low frequency trend in Figure 13 is due to altitude variation. The interband delay being a differential date phenomenon, its susceptibility to a sine perturbation depends on the time delay itself. Figure 14 shows the impact of this attitude perturbation for different band delays. Note that due to the linear susceptibility, the peak to peak perturbation relative to the delay is 0.15% for all couples. We have explained in section 3.3 that the yaw steering guidance induces a yaw drift. The consequence is that the west and east edges of the swath do not share the same ground velocity. The difference is proportional to the yaw angular speed, which depends quasi linearly on latitude ( Figure 5). Figure 15 shows the empirical ground speeds variations computed for center of detectors D01 and D12. In order to cancel altitude modulation, we compare the ground speed to the nadir one. The expected opposite linear trends are visible despite the noise coming from vibrations. Figure 15: Variation of Vground / Vground_nadir with latitude for center of extreme detectors D01 (blue) and D12 (orange). The opposite linear trends are due to the yaw drift which depends on latitude.

TIME DELAYS ESTIMATION FROM ORTHORECTIFIED L1C PRODUCTS
Now that delays are characterized in detail, we see that ESA interband table (Table 2), which is up to now the only way to retrieve these delays, may not be accurate enough to some applications due to the fact that fixed values are given independently of the detector number or time (Yurovskaya, 2019). We are hereafter proposing an interband delay retrieval algorithm for L1C Sentinel-2 products, which is the first level of product distributed to end-users. The section describes the algorithm and then presents the error assessment.

L1C products
L1B products are tiled orthorectified geolocated images (Martimort, 2007). Contrary to L1B products, there is no rigorous geometric model available so that methodology presented section 4 is not applicable. L1C tile size is 110x110km. Up to 5 detectors can lie on one tile. Thanks to the detector mask included in the metadata, users can identify each detector's footprint in the tile. This mask allows to perform pixel-specific delay estimation.

Proposed method
We have seen in section 3 that the delays can be divided in ACT and ALT components, the first being static and the other dynamic. None of these two components prevails on the other one, therefore we wish to correct both.
The first level of correction is to assess static ACT delays. Then a second level of correction is done taking into account the satellite altitude at the tile's date. Indeed, this is possible thanks to the presence of ITRF position ephemeris in metadata of L1C products. The whole processing chain is illustrated in Figure 17.
First, the static ACT rigorous estimation is done on a single L1B reference product. This estimation is done following section 4 methodology for a single date and the chosen band couples. The result ∆ OP is stored in a lookup table that is meant to be distributed to end-users. The result is identical to Figure 11 dots.
Then, we wish to compute Hsat and values associated to that tile. The ITRF ephemeris of the datastrip embeded in the L1C metadata allow to compute the nadir ground track for each sample on ellipsoid (or at any wished altitude). This can be done thanks to ITRF to WGS84 transform. The satellite altitude ephemeris Hsat[k] is one of the outputs of a WGS84 transform, and [k] ephemeris can be estimated by differentiation of nadir position ephemeris. The ephemeris sampling rate is 1Hz. Finding the best ephemeris without interpolation is accurate enough for our application because the altitude variation is very low frequency (see Figure 13).

Figure 16
: Nadir sample selection associated to a L1C tile. Now the problem resumes in finding the right ephemeris which corresponds to the tile acquisition date. One can see on Figure 16 that the nearest sample to the tile center location yields a fair estimate of the tile date thanks to the yaw steering condition. Retrieved accuracy on Hsat and is far from enough. Indeed, the detectors ALT footprint being 34km long with a mean ground velocity of 6700m/s, the expected maximum error is 6s. Referring to Figure 13, the delay error is much less than the vibration error. One may choose a nadir sample for each pixel of the georeferenced image. Nevertheless, computing one value of Hsat and for the entire tile still yields negligible errors.
Computed an inter-band delay for the tile is then straightforward.
Since delays are linearly dependant to / ratio, the Hsat_ref and _ OP values are retrieved from the L1B reference product and attached to the delays lookup table. The estimated time delay is then: Figure 17: L1C interband delay estimation processing chain.

Performance assessment
The performance assessment has been performed on the same L1B image dataset presented in section 4.2 and on the same ground points. As for L1C products, the satellite position ephemeris are embeded in the datastrip metadata, which allows us to use the same methodology detailed in section 5.2 for Hsat and determination. The reference product is located 67°W,19°N. This choice of reference will be explained later. Contrary to Figure 11 the lookup table is built with only the center pixel of the detectors, which makes a total of 12 points per band couple. Control points are used in order to assess the performances for the whole field of view. These control points are the one measured in section 4. They are located at the two detector's edges plus detector's center. Following Figure 17, for each ground point we compute Hsat and through the satellite position ephemeris. Then the lookup table is used with the corresponding detector in order to compute the estimation of ∆ . The value is compared to the one found in section 4 with the rigorous model. Results are presented Figure 18 with B2-B4 couple. Since the results are expressed relatively to the mean B2-B4 delay, it allows readers to extrapolate the results to other band couples. The relative errors are sorted by detector, column number in the detector, and latitude which allows to highlight the dependencies.
Compared to residuals obtained with ESA handbook table, we clearly see that the distortion pattern has vanished. Without surprise, biases are present at the extremities of detectors (light blue in Figure 18), since only the center of the detectors has been calibrated in the reference table some residual distortion remains. Our choice to calibrate only the centers of detector is motivated by the fact that it is laborious to evaluate the L1C pixel position in the original focal plane detector from the L1C detector mask. The loss of accuracy due to this choice is maximum 0.5% for edge of field detectors and negligible for center of field.
The linear dependency with respect to satellite altitude has also been very well corrected, as Figure 19 suggests. We note in Figure 18 and Figure 20 that the latitude dependence due to yaw steering has not been corrected and is one of the main limitation. Indeed, the proposed method do not correct for this kind of dependence. It could be done with a complementary lookup table or an empirical polynomial law. The bias is null for the latitude of the reference product, i.e +19°. That is why we chose a reference product that is located at mid-latitude. Being aware of the biases induced by latitude, one may choose a reference product nearer to user's needs. This latitude dependence yields a maximum error of 0.5% at high latitudes.
Finally, the performance of the proposed method depends both on the position in the field of view and on the latitude difference with the reference product. The maximum error is -1% for right edge of D01 detector at 80° latitudes. The overall rms error is 0.28% considering 3 points per detector and 0.17% considering the center of detectors. Using ESA handbook table yields a max error of 2,5% and a rms error of 1%.   detectors. These errors are due to the fact yaw steering speed has not been taken into account.

PERSPECTIVE ON A REAL-WORLD APPLICATION
The question might arise to what extent an incorrect time-delay plays a role in real-world applications. To investigate this let us consider the case of bathymetry estimation in the coastal zone using wave kinematics that deploys small time differences between satellite images (Poupardin et al., 2016 or image bands , Bergsma et al., 2021, de Michele et al., 2021. In these cases, the linear dispersion relation for free surface waves is used to inverse a local depth from ocean wave properties. In which h is the water depth, c is the wave celerity, k represents the wavenumber and g the gravitational acceleration. We impose an arbitrary wave period of 10 seconds, but for the analysis, in a dimensionless form, this is not important. We consider the minimum and maximum time lags between B2 and B4 found in section; respectively 0.986 s and 1.025 s while we use the ESA handbook table value for the time delay of 1.005 seconds as our base reference. In Equation 13, the time delay plays a role in the celerity (c). Except for the depth error in percentage, we express our analysis in a dimensionless form, on the y-axis hk, and on the x-axis Γ which is c²k / g. Figure 21: Error extremes for a depth inversion using the linear dispersion relation for free surface waves in relation to an inaccurate time-lag.
From Figure 21 one can deduct the maximum error related to the case presented in section 4.3. By only considering an inaccurate time-lag implementation, the error one makes for inversing the water depth is around 10 % for typical validity range of the linear dispersion relation (13), for 0.3 < Γ < 0.9. 10% is a significant error that can partially be accounted for by the proposed time-lag correction presented here, that typically reduces this error to a few percent.

CONCLUSION
In this study we have proposed analytic solutions to the interband pushbroom time delays with a detailed analysis of the contributions. Thanks to accurate estimations with rigorous geometric models, we have demonstrated the consistency of the analytic analysis, and put in evidence the true time lag variations. The optical distortion and the satellite altitude variations are the two main causes of the delay's dispersion. Finally, we have proposed an original way to retrieve inter-band L1C delays without rigorous geometric model. The performance assessment shows that we are able to estimate delays with 0.28% of rms error and 1% max error, which is at least 3 times better than the use of ESA table.
Enhancements are still possible: distortion residuals can be cancelled if one is able, thanks to the detector's mask, to estimate L1C pixel's position w.r.t the detector. Latitude dependency because of yaw steering can be mitigated through the use of a polynomial correction. We believe that a 0.2% peak to peak accuracy is reachable. This limit is the one due to LOS vibrations which cannot be easily mitigated.