Decennal geomorphic transport from archived time series digital elevation models: a cookbook for tropical and alpine environments.

—On the seasonal time scale, for accessible locations and when manpower is available, direct observations and ﬁeld survey are the most useful and standard approaches. However very limited studies have been conducted on direct observation at the decennial to century time-scale due to observational constrains. Here, we present an open and reproducible pipeline based on historical aerial images (up to 70 yrs time span) that includes sensor calibration, dense matching and elevation reconstruction over two areas of interest that represent pristine examples for tropical and alpine environments. The Remparts Canyon and Langevin River in Reunion Island, and the Bossons glacier in the French Alps share a limited accessibility (in time and space) that can be overcome only from remote-sensing. We reach a metric to sub-metric resolution close to the nominal images spatial sampling. This provides elevation time series with a better resolution to most recent satellite images such as Pleiades over decennial time period.

Abstract-On the seasonal time scale, for accessible locations and when manpower is available, direct observations and field survey are the most useful and standard approaches. However very limited studies have been conducted on direct observation at the decennial to century time-scale due to observational constrains. Here, we present an open and reproducible pipeline based on historical aerial images (up to 70 yrs time span) that includes sensor calibration, dense matching and elevation reconstruction over two areas of interest that represent pristine examples for tropical and alpine environments. The Remparts Canyon and Langevin River in Reunion Island, and the Bossons glacier in the French Alps share a limited accessibility (in time and space) that can be overcome only from remote-sensing. We reach a metric to sub-metric resolution close to the nominal images spatial sampling. This provides elevation time series with a better resolution to most recent satellite images such as Pleiades over decennial time period.

I. MOTIVATIONS
M OUNTAINOUS landscapes under tropical and alpine environments share a number of characteristics. Both environments show steep slopes, and their evolution are mainly dictated by climatic forcing (i.e., temperature, precipitation and extreme climatic events) which influences underlying mechanisms of geomorphic transport (e.g., soil formation, river dynamics, slope stability and mass wasting). Both experience important gravity-driven sediment transport (mass wasting, mass transport via glacier. . . ) that occurs from seasonal to decennial time span. Understanding how these environments evolve over decades is therefore essential to better anticipate developments in the 21st century in a context of global climate change.
As a matter of fact, the Intergovernmental Panel on Climate Change stated in 2013 that "more frequent and/or intense heavy rainfall events" were to be expected as a result of climate change "over most of the mid-latitude land masses and over wet tropical regions" [1]. In the one hand, hurricanes and storms generate destructive floods and mass wasting in tropical regions [2], [3], [4], leading to environmental degradation [5], [6]. In the other hand, in the alpine region, permafrost degradation and subsequent mass wasting occurrence clearly exemplify that phenomenon [7]. The consequences of these catastrophic events are amplified by other anthropogenic perturbations, such as deforestation or soil degradation. It re-Manuscript first submission on Feb 12, 2021 at IEEE GRSM. Revision awaiting editor decision. Corresponding author: A. Lucas (email: lu-cas@ipgp.fr). sults in biodiversity loss and ecosystem dysfunction [8]. In addition The remobilization of landslide deposits by heavy rains may have a devastating impact on the infrastructure and the housing, particularly in densely populated islands. For example, from 2004 to 2013, mass wasting processes were responsible for about 4,600 fatalities, and about US$15 million are annually spent for refurbishing damaged roads in the Caribbean [2].
Similarly, since mountain glaciers are highly sensitive to global and regional climate change [9], [10], sudden events like ice avalanches, outbursts of glacier-dammed lakes or floods on moraines can lead to large sediment loads in mountain rivers and to hazardous events with many casualties and extensive destruction. In addition, with the degradation of the permafrost, mountainous areas are expected to experience more mass wasting events in the next decades.
Therefore, as a result of climate change, tropical and alpine regions are faced with major economic and social challenges [11], [12]. Furthermore, gravity-driven processes (i.e., mass wasting and glacier) are some of the most efficient events that shape the surface of the earth. With underlying mechanisms that operate from a few seconds to 1000's of years and from cm 3 to 10 6 m 3 , they represent a major agent of erosion and sediment supply to rivers [13].
In light of these considerations, an assessment of landscape degradation over several decades, at the watershed scale and with a metric resolution is necessary. Measurements of changes in topography have long been used to quantify landscape evolution, and the emergence of satellites and remote sensing improved the spatial resolution and facilitated the studies on small timescale. Many studies have used a 4D approach to access the spatio-temporal dynamics of geomorphological processes at the decade timescale. Some of them used the association of various data sources for deriving topographic evolution. For example, by mixing the use of in-situ topographic measurements with (i) photogrammetry derived from satellite images [2], or with (ii) the extraction of base level from orthoimages obtained by satellites in combination with one digital surface model (DSM) of reference [14], [15]. In the one hand, despite a metric spatial resolution for the former method, such approach is limited in both temporal and spatial coverage due to limitations on field accessibility and/or the requisite workforces for in-situ measurements. On the other hand, the second approach is also limited, particularly by the resolution of the DSM of reference (i.e, ∼25 m in the case of typical national services such as the USGS in the US, the IGN in France or SwissTopo in Switzerland).
Other works have mixed DSMs derived from historical images with DSMs derived from recent satellites images. However, differences in the types of source generally lead to a significant thresholds in detection of change (typically >15 m) [16].
Additional studies have used times series from satellite images with stereo capabilities [17], [14], [18], [19]. Such approach offers high spatial resolution, and has been welldeveloped and intensively used since the 2000's as it is well supported by commercial softwares. Similarly, recent sensor developments on LiDAR (Light Detection And Ranging) and UAV (Unmanned Aerial Vehicle) make very precise measurements even for the ground surface below the canopy in the case of LiDAR. However, both methods are only limited to the post 2000's period.
Finally, one can take advantage of archived images that offer a great historical depth as well as a sub-metric spatial sampling. Nonetheless, the use of original photographs implies intense manual edits, and until now, most studies have dealt with a narrow spatial cover on specific objects such a individual glacier or individual landslide, that requires small enough number of images to be manually processed [20], [21].
Although these different methods are all applicable for geomorphic processes analysis, the question of reconciling large spatial coverage with a multi-decades time coverage still persists. To answer this problem, the most suitable approach would be to compare time series DSMs derived from archived aerial images, that allow to assess topographic information over several decades (pre-2000) with large spatial coverage (>10 km 2 ) and high resolution (sub-metric).
Nonetheless, the main difficulty in such an approach lies in the fact that historical aerial photos have been acquired by film cameras and that even if they have been subsequently scanned (with risks of bad conservation over the years) they were not acquired digitally, leading to flaws and subsequent difficulties in deriving accurate time series DSMs from photogrammetry algorithms. Moreover, depending on the time and space scale, the quantification of geomorphic processes by photogrammetry requires a large number of images. In spite of the existence of software or photogrammetry programs (free or proprietary), that are very powerful (and easy to use for recent images), the use of scanned archived photos requires a very important preparation to be able to use the algorithms of photogrammetry. In that regard, the use of a large number of images for detailed time series, large geographical areas and detailed geomorphic analysis, requires the automation of all the steps that can be automated, for the production of DSM from archived aerial images.
In this paper, after highlighting the limitation of current methodologies for DSMs production when using archived scanned images, we introduce a new workflow with automatic steps and we show its application and the results on the production of DSM for tropical and alpine contexts, namely Remparts and Langevin rivers (Reunion Island) and Les Bossons (french Alps).

II. FROM SCANNED AERIAL IMAGES TO DIGITAL SURFACE MODEL, AND THE ASSOCIATED PITFALLS
Automatic processing archived images for photogrammetry is a difficult task because numerous steps are hampered by flaws that can make automation very difficult and/or propagate errors throughout the entire workflow.
To synthesize, digital photogrammetry applied to overlapping aerial scanned images is based on 4 main steps. The first one consists in recovering, from the raw images and from external sources (e.g. camera calibration report, for sensor and optical system), the technical information necessary to use the scanned photos as if they had been acquired with a digital camera (e.g. digital single-lens reflex camera (DSLR), UAV). The second step consists in extracting the relative orientation between all the images, using tie points (determined either manually or automatically, depending on the software considered). The third step consists in placing this resulting orientation in a georeferenced (i.e. absolute) coordinate system using ground control points (GCP) acquired either in the field or from a reference basemap. Finally, in the fourth step, pixel correlation is used to derive the final DSM for the area covered by the aerial photos. In detail, each of these four steps requires particular caution and are associated with limitations that we highlight below (by being as agnostic as possible with software).

A.
Step #1 and the potential issues regarding the assessment of autocollimation Before being accessible as downloadable images from online catalogs, the scanned images were acquired from a socalled "pinhole camera", and after being physically archived, they were digitized using a scanner. Therefore, an important task is to process these pictures to obtain images as if they had originally been taken by a digital camera, and to obtain the related metadata (i.e., this what is called pre-processing in the following sections). First, it is necessary to find the position (in the scanned photo) of the fiducial marks in order to calculate the exact internal orientation of the pixel in the images. Even if this task is easy for the human eye, and in order to make it automatic to process more than tens of photos, it is necessary to take into account several aspects. For example, the pattern as well as the number of these marks can change from one camera manufacturer to another, and when it comes to multidecennial archives, many patterns must be taken into account. In addition, as shown in figure 1, photos can be in grayscale or color ( figure 1-a,b). In some cases, fiducial marks may even be barely detectable when covering low contrast areas (e.g., figure 1-c). Others can be small ( figure 1-d,e). Finally, fiducial marks may be either cut or missing due to misalignment of film in the scanner during the manual digitization procedure (e.g., figure 1-f).
Once being identified and located in each pictures, these marks are used to re-project and crop the images so that they are mapped as if they had been acquired digitally in the first place. For this remapping, an affine transformation is generally considered to overcome potential misalignment due to the manual scanning procedure, After this remapping it is necessary to define the intrinsic geometry of the sensor (height and base of the field of view pyramid that corresponds to the photos) using the specifications of the camera such as (i) the focal length F , defined as F = f i|j × p i|j , with f i|j the focal length in the i, j axis, and p i|j the pixel size in the i, j axis, (ii) the main point ([i 0 , j 0 ]) called boresight, and (iii) s the skew coefficient of the axis that causes shear distortion and aperture or sensor/pixel size.These three parameters are used to derive the intrinsic matrix (i.e., internal calibration) of the camera K as: These parameters and specifications of the camera are often provided with a calibration certificate but they can also be reestimated, although it usually requires additional calculations for a good convergence.
This procedure of cropping/re-mapping the pictures and of deriving the internal calibration of the camera is critical. Each photos of the same survey being taken with the same camera, the the matrix K is constant for all the photos of this survey, and to be relevant, pictures must be accurately cropped along the fiducial marks.
Depending on the software package used for the autocollimation, the steps described above can be purely manual or partially automatic. In the case of a manual procedure over a large data sets, in addition of being time consuming, it may result in pixel misalignment between images due to inaccurate pointing of fiducials marks or inaccurate cropping, etc. Such inaccuracy has a direct impact on the quality of subsequent operations, and thus on the resulting products. In the case of semi-automatic procedure, when it exists (i.e., at authors' knowledge, no purely automatic procedure exist nowadays), we still observed (in the MicMac suite, e.g., [22], [23]) up to 10% of misalignment due to erroneous fiducial marks detection over more than 100's of images of the same survey.
This first step of auto-collimation results in the internal calibration of the camera, that will be used to estimate the relative orientation of the camera that took the pictures under consideration.

B.
Step #2, the relative orientation This step called relative orientation consists in inferring the relative position of the sensor for every images, accounting for (i) the homologous points (tie points) between images and (ii) the internal calibration of the camera obtained from the previous step. Depending on the type of implementation (that depends on the software package), this step may also be used to derive internal parameters, such as the optical distortion from the optical system on board.
This procedure of relative orientation is automatically performed in an iterative way and starts first with the identification of tie points along each possible image pairs. However, this operation using a feature detection algorithm for image analysis (i.e., SIFT in the workflow presented here) is time demanding as it goes in n * (n − 1)/2, n being the number of images to process. To facilitate and accelerate the procedure, one can provide a reference list of image pairs. Once enough tie points are detected for each overlapping pairs, a second operation consists in retrieving the relative orientation of each camera. This is done from inferring the sensor matrix defined as: where R and t are the extrinsic rotation and translation, respectively. This step is also called "external calibration" and sometimes integrates the estimate of the radial distortion of the camera optical lens system (for some software it is performed in an intermediate step). In order to take into account the optical distorsion in the position of the pixel on the image, it is necessary to estimate the radial distortion coefficients k n , used to calculate the undistorted coordinates idx corr (i.e., index in the vector image) as follows : where, r is the conic, defined as r = i 2 + j 2 . In practice, for long-focus lenses, like the ones used in aerial surveys, two coefficients of the radial distortion should be sufficient.
The result of this step is the relative arrangement of each image. Yet this relative orientation does not fully represents the reality of the scene one seek to reconstruct since it is not georeferenced or even scaled, which is done in the next step.
However, it is worth noting that in order to facilitate the procedure of tie point recognition, one can directly provide the approximate locations of the images. Indeed, depending on the national service responsible for aerial survey, flight path are often available and so provide embedded GPS-like positions. Nonetheless, for campaigns prior to the introduction of GPS (early 1980's) we observed errors from 10 to ∼100's of meter in horizontal locations between position approximated from the flight paths and the final DSM (accurately georeferenced). Yet, elevation, obtained from the plane altimeter is much more accurate even for old campaigns. As a consequence, the heterogeneity in the accuracy of the location of the images must be taken into consideration when using flight path. Finally, in this case where the locations of the images are available (i.e. a priori), this step of external calibration leads to an approximation of the absolute orientation, which is described in the next section.

C.
Step #3, absolute orientation and the crucial aspect of the Ground Control Points After computing the relative orientation for every image, another procedure is required in order to accurately place this orientation into a georeferenced coordinate system. This is done by retrieving the absolute orientation of the cameras (i.e., inference of R, t, w) using external geographical information, called reference points or Ground Control Point (GCP). Indeed, external control are required to scale up the sensor coordinates (row and columns of the pixel in the image) towards the world coordinates following: where [i, j] are the pixels in the image, [x, y, z] the points of reference in the world and w a scale factor.
Although for several decades, modern sensors have been equipped with integrated GPS providing geographic information during the acquisition and facilitating the transformation of pixel coordinates to world coordinates, archived images do not provide such metadata. Therefore, the step described here require GCPs that provide the exact geolocation M (x, y, z), of a particular site/object on the images. GCPs can be of various nature and generally require manual recognition or manipulation on the photos under investigation. For example, GCPs can be GPS measurements from field surveys, well-identified and recognized features on orthoimages or even benchmark points of well-known elevation selected on basemaps. In the literature, most studies acknowledge the usage of GCPs [21], [20], yet with the aim of studying DSM time series over several decades, finding GCPs over such periods is in fact difficult. The general approach is to manually recognize features that exist on the recent photos and on the old photos under investigation. Urban constructions have been widely used as GCPs since it is easily recognizable on pictures, but even for these features it may still be challenging to recognize them in old pictures, especially when urban areas have undergone significant development during the period covered by the studied images (extension of houses, enlargement of road junctions, etc.), see Fig. 2).
Of course, this is also true for features in landscapes with high rates of change (e.g., floods, landslides, rock avalanches, overflows), which is often the case for relief of steep slopes. Therefore, finding (on photos) and extracting GCPs that exist now and that existed before is not always an easy task and is a labor-intensive task, even for few images. For decadal time series, it usually results in a small number of GCPs (a few tens of points), sparsely scattered over the entire canvas or grouped in one place ( Figure 3). Indeed, there are two main reasons why GCPs cannot be recognized in campaigns of different ages. Shadows are one of the problems, while some GCPs may simply not exist on the photo of a different date. As a result, some areas are often over-weighted in the selection of GCPs (see the probability distribution in figure 3). For example, in Reunion Island (one of the two reference areas of this study and described later), the urban areas in which GCP recognition is easiest, are distributed along the volcano flank, resulting in a more or less homogeneous distribution of GCPs as a function of altitude (Fig. 3). On the contrary, in the Chamonix valley (the other area of interest described later), the agglomeration is localized in the valley, which gives a heterogeneous distribution of GCPs as a function of altitude. As exemplified, the difficult selection of GCPs, as it is a manual procedure (for all software), often introduces a strong bias in the altitude constraints and often leads to a number of GCPs lower than the actual potential number that could be used.
Finally, the relevance of the chosen GCPs is not only related to their use for the absolute orientation but also to the fact that they are used in a subsequent procedure to optimize the camera calibration (and thus the absolute orientation). Indeed, after the absolute orientation calculation, the resulting point cloud (which represents the studied scene) is used as a reference canvas from which, using the GCPs, the camera parameters are refined. This procedure is called the bundle adjustment. This adjustment may concern the relative locations of the cameras, the orientations of the cameras or even the coefficients of radial distortion (only when numerous GCPs are available). This optimization therefore depends heavily on the quality of the GCPs, which may ultimately, for the reasons mentioned above, be insufficient to obtain an accurate aerotriangulation of the whole canvas.

D.
Step #4 and #5: 3D points cloud and digital surface model reconstruction The final step consists in computing the depth map from ray intersections by using dense correlation over the images. This is done using image correlation from moving windows that are easy to implement and well suited for long range acquisition (i.e., small B/H ratio). However, archived images often show low radiometric dynamics that can lead to inaccurate DSM construction due to the saturation of very bright areas (i.e. glaciers) and/or to low brightness and/or low contrast in shaded areas (i.e. valley bottom). The overall workflow is summarized in Figure 4.
E. Potential issues regarding the quality assessment of the resulting DSM using GCPs Considering the criticality of each of the four main steps in the process of construction of DSMs from scanned aerial images, an accurate quality assessment of the result is required. The simplest way of estimating the accuracy of the resulting DSM would be to use the GCPs used to compute the absolute orientation and compare their position (x,y,z) to their corresponding location in the modeled DSM. In other words the accuracy of the resulting DSM could be inferred by the standard deviation of errors at the GCPs locations. However, if such a comparison is necessary, it is not sufficient and might sometimes be misleading instead. Indeed, although the step of absolute orientation (i. e. bundle adjustment) seeks to minimize the Euclidian distance between ray intersections and the location of the GCPs, too few and too sparse GCPs (over the landscape of interest) can lead (depending on the free parameters that are considered during this step) to a global distortion of the resulting DSM, while close to the GCPs [24]. Therefore, in the case of two campaigns that share little to no GCPs, such potential and global distortion will lead to a large difference between DSMs while having a small standard deviation of errors at their respective GCPs locations. The direct consequence is that such a difference (in altitude) between two DSMs will account for both the error in the DSM construction and the actual changes in elevation due to the real geomorphological processes that one seek to quantify. In addition, attempts in correcting a global distortion due to sparse GCPs (or to an erroneous estimates of the radial distortion associated with quasi-parallel sights [25]) using polynomial regression, should also lead to inaccurate final DSMs and then hinder for potential real changes in elevation. A strong mismatch between the GCPs and their equivalent on the final DSM clearly show the poor quality of the final product, but a good correspondence does not necessarily implies an accurate DSM, espacially with sparse GCPs. Other tests such as the the variation in elevation over areas known to be stable over the decade time scale studied need to be evaluated.

III. COOKBOOK FOR A ROBUST PHOTOGRAMETRIC WORKFLOW FROM HISTORICAL TIME SERIES
In order to automatically produce time series of DSM, from large volumes of archived images and within the previously identified limits, we propose a complete workflow with new approaches on pre-processing and canvas co-registration steps. In this work, our efforts are mainly focused on sensor calibration since it is the most critical and difficult aspect when processing scanned images.
A. Pre-processing As described above, a crucial step is to transform the scanned films into digital type images using the reference marks. While many software packages offer no solution, some others (both commercial and open-source) have rudimentary tools that essentially consist of locating the marks by clicking on the images. Depending on the software, either the manual Step-2-

Point cloud construction
and/or camera position pointing must be repeated on each image, or the spot automatically propagate across the data set assuming a constant or at least very close location for all images. By applying the latter method with MicMac [23], we have observed more than 10% of errors on more than 100 images, which then have to be corrected manually. As long as the fiducial marks exist in all scanned images, the next step of co-registering the reference marks can be achieved with a sub-pixel accuracy. In addition, the manual procedure of scanning archived images often leads to missing and/or partially cropped fiducial marks, that prevent the manual or existing automatic identification process. To overcome these limitations, and to allow processing hundreds of images from various campaigns with heterogeneous fiducial mark patterns (see Fig. 1), and potential missing marks, we introduce the FiducialLib workflow as outlined in the figure 5.
This workflow consist of four procedures. 1-The Fiducial marks pattern extraction is a manual step that defines the template images corresponding to the types of fiducial marks present in the images of the different campaigns under consideration. 2-Manual creation of a catalog that describes the number of fiducial marks in the images of the different campaigns (i.e 4 or 8 marks). 3-Automatic search all the images for the fiducial marks by pattern recognition based on a simple correlation operation: where S is the score of the pattern recognition step from the fiducial mark pattern template F into the image I, and where prime function is defined as: where w and h being the width and the height of the pattern template. The location of the fiducial mark corresponds to the best match, and is being obtained as the global maximum of S.
In the case of a missing fiducial marks (due to the cropping of the image during the scan procedure), a virtual frame is modeled from the other identified fiducial marks, so the location of the missing marks can be retrieved, in order to reshape/remap the image accordingly. Assuming a 2D homothetical deformation, a minimum of three fiducial marks are required.
From the whole data used to develop the workflow presented in this study, we empirically found that keeping max(S) > 0.5 allows us to retrieve the different pattern templates even when they are located over an area of a very different brightness in the different images (e.g., forest vs. snow).

Fiducial mark location
Numerical-like image data base Finally the fourth step of the FiducialLib workflow is the remapping using a homothetic transformation. This operation takes pixels from one place in the image and locating them in another position in a new image resulting into a digital-like image, as it was taken with a digital camera in the first place. The whole FiducialLib is provides in algorithm 1.

Image remapping from affine transformation
An example of the use of the FiducialLib is given in Figure  6. The main advantages of our approach is that it is an agnostic method with regards to the fiducial mark patterns and needs very little manual edits. After the manual definition of the template (as many as the number of campaigns required) Fidu-cialLib automatically process the images. More importantly, over the +100's images processed for this study (i) we do not observe any rejection over various campaigns, even when one to two fiducial marks are missing, (ii) we achieve a sub-pixel accuracy for the locations of the fiducial marks, and (iii) the missing mark problem is overcome as long as two orthogonal axis can be identified.

B. Canvas co-registration
In order to produce accurate time series of DSMs, all considered canvas from the different campaigns need to be co-registered onto a common georeferenced frame. Yet, the second major limitation of photogrammetry using historical images to produce time series DSM, is associated to what we name "the GCPs paradigm". While required for bundle adjustment and georeferencing of each canvas, manually chosen GCPs do not necessarily guarantee an accurate co-registration between canvas from different campaigns, essentially due to the limitation discussed in the previous section. The ideal approach to maximize the accuracy of the co-registration of the canvas would consist of selecting a large number of GCPs well-distributed over the area under investigation and which exist through all the canvas. Unfortunately as discussed in section 2, this is barely possible in real situations. Therefore, another approach consists of selecting different GCPs only on the canvas of reference C 0 . The two main advantages of using the most recent canvas as the canvas of reference (to register the others) is that (i) it allows to maximize the chances of selecting and recognizing numerous GCPs from other reference sources such as orthoimage, maps or GPS measurements on the field, and that, (ii) if the most recent canvas is acquired from digital camera, it helps optimizing the sensor calibration convergence. Yet, although selecting GCPs is essential for the absolute orientation and the bundle adjustment of each canvas of a time series, their choice, their nature, and their use in the alignment of the canvas are critical. Consequently we propose to avoid the use of GCPs to align the different campaigns (i.e., canvas) with each other, by using, instead, a feature detection algorithm such as scale-invariant feature transform (SIFT).
Indeed, even in a changing landscape, many areas or objects remain stable over time, and the number of stable features in an image increase with its coverage area. For instance, human recognizable objects such as fixed roads, old buildings, and some crests and hills do not necessarily evolved over the time frame considered (e.g., 5 to 10 years which is the typical time lapse between aerial survey). But more importantly, feature detection algorithm does not deal with real objects but with image derived features. It extracts its own feature description and use it to identify the same feature in another image, thus overcoming eye recognition problems and thus ensuring a better spatial distribution and a greater number of these tie points between each campaign and the reference one C 0 .
By doing so, we obtain a sufficient number of homologous points features in a common scene, typically between 100's and 1000's when images overlap with more than 50% in coverage.
This step is done by applying a pre-processing difference of Gaussian filter on the images, that allows to catch some point even in either saturated (i.e, glacier/snow cover) of under shadow (steep cliffs) areas, which are common on our region of interests. In addition, in our approach there is no need for all images of a given canvas C n , to be link/associated to the canvas of reference C 0 as long as each image of canvas C n is associated with at least one other image of C n . The geographic information from the reference C 0 is hence propagated through the entire series of C n , backward in time from one canvas to another, with C 0 the most recent canvas.
Moreover, we take advantage of this approach for optimizing the co-registraion of the different canvas from the different campaigns. Indeed, processing archive aerial images usually lead to bad convergence during camera calibration, that can be improved by setting the sensor properties (focal length and boresight, when initializing the calibration). However, since the different canvas C n are all related by new GCPs (tie point between campaigns), we can propagate the orientation information from C 0 through C n . As the canvas of reference needs to be the most recent campaign (see above) the orientation is propagated backward in time, from one canvas to another with no re-evaluation of the reference canvas (i.e, C 0 to C 1 , C 1 to C 2 etc. or C 0 to C 1 , C 0 to C 2 etc). In practice, this is done by performing simultaneously both, the external orientation and the bundle adjustment with 100's to 1000's tie points for each image. As far as the spatial distribution of the tie points (between each campaign and its reference) is not statistically biased, then the process converges rapidly and leads to co-registered canvas after the step 1 to 3 as described in section 2. In a similar manner, we can backpropagate the aero-triangulation to older canvas. In the case Algorithm 1 FiducialLib's algorithm Input: Fiducial Mark F + sensor canvas SC + I 1 · · · In scanned images Output: J 1 · · · Jn digitized-like images return J 35: end function of a too important age difference between two campaigns, one can chose a given C n becoming a new reference for a "very" old canvas. We call this approach the "backward time nearest neighbor propagation" (BTN2P). Finally, since it relies on a well established bundle adjustment algorithm, the accuracy of the canvas co-alignment is easy to assess.
It is worth noting that, an multi-block (or multi-canvas) approach has been proposed for painting restoration using multiple cameras [26]. However our methodology varies in many regards. First, the scene depth for landscapes are changing overtime, while in the paintings case, changes are in mainly in the radiometry (i.e., the restoration does not affect the depth of the scene). Second, historical surveys come with a badly constrained orientation or no flight path, while in the case of painting restoration, the location of the cameras are accurately known and fixed. Finally, and being the most critical point, the scanned sensors are badly calibrated as discussed previously. Because of these considerations, in our method, both internal and external calibrations are done by the "backward time nearest neighbor propagation", as opposed to an evaluation by iterative addition of canvas [26] that would lead to non-convergence in the case of archived scanned images.

C. Overview on the workflow
Here we summarize our new approach for the production of time series DSM from archived aerial images: -I-Download the scanned archived aerial images from online catalogs or manually scanned achieved films, of each campaigns of the studies times series -II-Pre-process all the images from all the campaign with FiducilLib (algorithm 1) in order to transform the scanned images into images taken by digital cameras.
-III-Designation of the campaign of reference and selection of the GCPs for the campaign of reference. In our strategy the archived images the closest in time to the dataset (e.g. maps, GPS measurements) used to select the GCPs, is the campaign of reference. Selection of high quality GCPs on the images of the campaign of reference from basemaps of national surveys/institutions (e.g. IGN for France, SwissTopo for Switzerland, USGS for USA...), orthoimages or GPS measurements.
-IV-Process the campaign of reference and evaluate the relative and absolute orientation of its cameras (Step 2 and 3 of the classic photogrametry workflow). The result is the canvas of reference, named C 0 (function ImgTo3DCloud(C 0 ) in algorithm 2).
-V-Process the other campaigns using the "backward time nearest neighbor propagation" approach. The tie points between the considered canvas C n with C 0 are used for both evaluating the internal and absolute orientation of the cameras (function ImgFromC0To3Dcloud(C 0 , C m ) in algorithm 2). Note that the evaluation can be performed from an orthoimage as far as a common tie points can be found with C m . -VI-Dense pixel matching resulting in 3D cloud for each resulting canvas C n using shape from motion algorithm upon user choice.
-VIbis [optional]-While this does not concern the data discussed here, in case of poor tie point matching between canvas, and/or in case of unbalanced spatial dispersion of the tie points between canvas (as it may happen when radiometry varies between canvas and or if the time span is too large in regards to the actual changes in the landscape), hence a 3D homothetic compensation can be done by assessing the Euclidean isometry that minimizes the point to point distance between the considered 3D cloud and a reference considering a residual threshold upon user choice (either the 3D cloud derived from C 0 , or from an external source, function CloudCorr(3Dpts 0 ,3Dpts m ,Threshold) in algorithm 2). The computing transformation matrix is applied to external orientation of the cameras, so both 3D cloud and derived orthoimages can be corrected with this Euclidean isometry in order to preserve the relative shape of the cloud -as opposed to polynomial correction -allowing mass balance calculation of time series.
We provide the whole workflow in algorithm 2.

A. geomoprhological settings
In order to demonstrate the benefits of the workflow presented in this study the analysis of the rate of geomorphic change, we applied our methodology to two steep-slopes environments from two different geological contexts. First we chose a region in the tropical environment of Reunion Island, in the Indian Ocean. The area of interest is groups the Remparts and Langevin rivers [2] that incised the flank of the Piton de la Fournaise volcano over the last 60 Ka [4], leading to a deep canyons between the remaining and stable parts of the volcano (called planezes). In these canyons, the difference in elevation between the planeze and the valley bottom ranges 1500 m to 90 m with valley slopes up to 60 degrees. The second area of reference is the Bossons glacier located in the French Alps, next the the city of Chamonix. In this region the geomorphological context is different from the Reunion Island since glacial valleys are shallower and wider than river canyons. However the glacial valley des Bossons presents a steep slopes of flow of 20-30 degrees between 1200 and 3000 m. These two regions have been chosen for their opposite climatic environment, for their steep but different kind of relief and for their strong elevation change between the different geomorphologic parts. This particular aspect carcaterizes, what is call in photogrammetry, the scene depth which is known to be a technical challenge.

B. Data sets
The archived aerial images used in this study were acquired by the Institut Géographique National (IGN), the French national geographic service surveys. This office have provided surveys over the french territories since the 30's and have acquired aerial image with a frequency of about 5 to 7 years. Depending on the period, the flight conditions (elevation) and the camera focus length (80-220 mm), images present a ground sampling ranging from 10's cm to 80 cm. In this study we chose to process a number of 4 campaigns, from 1949 to 1988 for the Bossons scene, and from 1978 to 2003 for the Remparts-Langevin scene. The number of images per campaign ranges from 6 to 21, for a total of 37 images for the Bossons scenes and 67 images for the Remparts-Langevin scenes. The reason why the Bossons scene has less images is that the topographical settings of the Chamonix Valley does not allow any flights at low ground elevation. Therefore the footprint sized of the images of the Bossons scenes are wider compared to the Remparts-Langevin scene (acquired at lower altitude) and a fewer number of images is needed to cover a comparable area. Both scenes cover an area of about 100 km 2 . Finally, the digitized historical images were downloaded from the IGN web service, the size of the images ranges from 7000x7000 to 12000x12000 pixels, and depending on the period, images are in grayscale before 1984 and for some campaigns in colorscale after 1984.
Note that such a data sets required a heavy computing powerhouse, such as high-end workstation or HPC facility due to the highly CPU usage as well as memory requirements for some specific steps, such as the search for tie points and the dense correlation process.

C. Results
The landscape dynamics of these two areas is out of the scope of this study where we focused on the methodology. We implemented the two algorithms in a mixture of python codes and MicMac/Apero commands but they can be applied to a various kind of software at the reader's discretion. Example results over the two areas, obtained from our new workflow are shown in Figure 7. Each cloud corresponds to a region of about 80-100 km 2 with a median geodetic distance of 0.8 meter between the nearest point. As a matter of fact, spatial sampling of the resulting 3D cloud is directly a function of both the image sampling and the image redundancy. Thus, the 3D cloud has a spatial sampling of a factor of 1 to 4 compared to the image. In our regions of interest, some areas are covered with up to 28 images. In such areas, the resulting 3D cloud has a spatial sampling similar to the image's. Consequently, in order to derive gridded DSM time series with homogeneous resolution, we decided to resample the 3D cloud at 1.25 m/post, which is superior to the median sampling, to guarantee independence between the nodes of the elevation grid. This aspect is critical when using the DSM to derive slopes, mass balance and transport fluxes etc.
For both of studied areas, we obtained residuals of pixel to sub-pixel RMS (0.3-1.2) accounting for the data set covering the 1939-2003 period. This assures a metric to sub-metric coregistration between canvas, considering a ground sampling of 15 to 80 cm.
As stated previously, the use of GCPs for the estimate of the accuracy of the resulting DSMs is somewhat misleading because such analysis does not account for any global distortion of the block. Nonetheless, the residuals at our GCPs fall in the metric order of magnitude in z, being taken on a DSM of reference with a ground sampling of 5 meters. Typical horizontal residuals are below one meter where the vertical residuals range from 10's cm (on flat or built terrain, such as houses) to 5 m (where there are picked close to steep slope badly resolved on reference DSMs, i.e., IGN's RGEAlti). Most importantly, the residual in between campaign at GCPs location fall below 1 m. But again, this good agreement at GCPs between campaign is not sufficient for assessing the ability of our resulting 3D cloud to be used for landscape evolution over some decades. Instead, we examine the quality of the results using stable geomorphic areas. Indeed, comparing the same stable areas between DSMs of different dates allows to estimate the errors due to the whole process of reconstruction of the DSMs, and allows to estimate the sensitivity of our approach and especially the use of automatic GCPs.
Although only a small part of any landscape does not change over time, stable areas increase as the period of time decreases. Consequently at a decennial time scale, it is possible to find, even in active regions, parts that are not affected by erosion or other geomorphic processes. In the case of the Remparts and Langevin rivers, the flank of the Piton de la Fournaise volcano (planezes) that are not affected by river incision can be considered as stable in a 50 to 100 yr timescale. In the Bossons Glacier region, the bare slopes of the glacial valley can also remain stable over a several decades. In order to evaluate the accuracy of the method, we performed the analysis of stable geomorphic zones in these two types of regions. The figure 8 shows the distribution of the differences in elevation in the stable areas (geomoprhic stable areas), between two campaigns relative to the elevation. Differences in elevation range from 5 cm/yr and 92cm/yr with little fluctuation of the median around zero in the whole range of elevation suggesting no particular trend in the variation of the difference of elevation with latitude, longitude or elevation. Consequently it shows that our approach results in an optimal calibration and orientation of the cameras all across the campaigns. As a matter of fact, we observed no problem of poorly established focal length, which is one of the most sensitive parameters during the internal calibration. Indeed, a bad estimation of this parameter leads to a erroneous depth of the global solution, hence to the slopes of the scene. Here, a good solution is obtained for both regions under investigation. Moreover, the absence of trend also shows the absence of relative tilt of the aerotriangulations of the different campaigns, whereas in our approach only one campaign is conditioned with GCPs. This highlights the relevance of our approach as being effective while being time saving. Finally, the figure 8 shows that in the one hand, for the bare slope of Chamonix Valley, the change in height of the modeled DSMs range between 1100 and 3000 m suggesting that our approach allows to detect geomorphic changes in the order of a few decimeters per year. On the other hand, in La Réunion, in the stable area that contains vegetation the the change in elevation range from 750 m to 2200 m leading to a limit of detection of our approach of 1 m/yr.
In addition, it is worth noting that, as we derived the surface model field from optical images, the resulting products combine bare rock terrain and vegetation cover, and hence in a time series, differences can show the vegetation evolution over time (e.g., agriculture fields, canopy changes etc.) in region geomorphologically stable (where bare rock terrain does not change).
To summarize, we have developed a method for generating DSM time series from archive scanned images. The results are a time series DSMs that are intrinsically co-aligned with each other with a relative precision of the order of one meter. Only one canvas (i.e., block) of the investigated time series needs the use of GCPs which dramatically reduce the associated flaws and limitations. Another advantage of the presented workflow is the reduced manual steps allowing the automatic processing of large data sets of images (i.e., 100's). Finally, the resulting precision allows to detect artefact due to the scanning of the archived images. Figure 9 shows an example obtained over a very stable part of Piton de la Fournaise lava flows between 1984 and 1978 where parallels bands of ± 0.5 m are clearly observable along image edges, indicating the direction of the scanning process.

V. DISCUSSION
This workflow is presented here as a cookbook for environments that show very steep slopes. The two different climatic environments used here to apply our method both share contrasted reliefs, but with different distributions of slopes because of different geological, morphodynamics and climatic contexts. In addition, urbanization and land use have evolved at different pace over the considered period. However, in both cases our workflow allows detecting and quantifying the morphodynamical and landscape evolution over decennial time period.
Notwithstanding, few remaining issues needs to be discussed and addressed. Some of them, related to the nature of the raw data, are unsolvable while others can be partially resolved.
Aerial image sensors operate in visible spectra and are sensitive to climatic conditions. Therefore, aerial photo often show partial to significant clouds coverage. This particular  The blue curve show the median yearly rate between two dates. The color map represents the probability distribution of the yearly rate for each elevation bin with non-black being 90% of the distribution. On the Bossons scene (top), we see the trees evolution around 1100-1600 m while most of the bare ground is stable above (excluding glaciers). The planeze ontop of Remparts-Langevin (bottom) are corresponding to 1750-2250 elevation range where minimum change rate is observed. The valley corresponds to the 500-1600 m elevation range. White areas on Bossons scene is due to smaller changes compared to Remparts-Langevin, when excluding glaciers and river respectively. point is critical for tropical areas as well as for high relief regions. This causes two unsolvable limitations: i) the clouds are not optically transparent, and ii) in case of partial cloud coverage, their shadows on the ground limit the DSM reconstruction. This is even more critical under windy conditions, when clouds are moving from one image to another during the flight acquisitions. At this stage, the best practice would suggest to get rid of such images.
Another issue that can be encountered is related to the topography of the investigated scenes since steep landscapes will lead to shadowed hillslopes. In the shadow, the radiometry range drops by a factor 2 or 3 compared to other areas and make the dense matching very difficult. Attempts to increase the brightness in such areas make possible to a better pixel matching. Nonetheless, topographic reconstruction of these areas remains less accurate than for the ones illuminated. We show, on the other side of the spectra that very bright areas such as glaciers, are well reconstructed because associated pixel are actually not saturated on the considered image data sets (e.g. Fig. 7-a).
Finally, since our workflow allows production of DSMs with metric resolution, and because biomass growth can be very quick in different environments, the method presented here can help assessing the canopy height evolution and hence biomass dynamics over decades in different climatic settings. This is a critical aspect of either land use science and/or ecology research, and we suggest that a workflow using archived aerial images can be used for other purposes and not only for quantifying geomorphological processes.

VI. CONCLUSION
A complete and self-consistent workflow for deriving time series DSM from scanned archived image is presented. The presented global and automatic method accounts for the various kinds of fiducial Mark patterns that exist over the various campaigns. It requires a minimal usage of GCPs that dramatically reduce the manual edit steps and associated limitations. The resulting block association and compensation allows to obtain well aligned canvas. As algorithms are provided, this is up to the reader to implement it in his/her software of choice. Consequently, the presented workflow can be use to obtain time series of DSMs over decennial time period with a sub-metric sampling. The additionally outcome of the methods is to provide accurate true orthoimage time series. This is especially of interest for time evolution assessment for landscape evolution experience climate forcing [27].
When applying the workflow to two highly steep terrains, surface process with annual rate change of the order of a meter are quantifiable. This allows to study landscape dynamics with sediment flux and mass balance assessment in the decennial time scale which is critical in geomorphology. The DSMs resulting from the method presented here are therefore useful for mass movements assessment, such as landslides, glaciers and the like. Finally, as we derive DSM, canopy height can be also derived, hence biomass evolution. Consequently, the resulting products of such a pipeline offer substantial new insights into the critical zone science.