Use of Earth observation for monitoring soil sealing trends in Flanders and Brussels between 1976 and 2013

The on-going growth of urban area in Flanders and in the Brussels Capital Region over the past decades has resulted in a highly sprawled urban tissue, consisting of large and smaller urban agglomerations, connected by a well-developed transportation network. The conversion of open land to urban area is accompanied by an increase in soil sealing, affecting the hydrological cycle and the urban climate. Despite a growing interest in monitoring the process of soil sealing in urban areas, to date no detailed information on the presence and evolution of sealed surfaces is available for Flanders. In this paper a linear regression unmixing approach is proposed to map and monitor changes of sealed surface cover at the regional scale, using medium as well as high resolution remote sensing data. Applied to Flanders and the Brussels Capital Region, a total sealed area of 2687 km² for 2013 is found, corresponding to an increase of 82% since 1976. Residential areas account for nearly half of the sealed area and show the largest increase in sealed surface cover over the past 37 years.


Introduction 1
Anyone advocating the compact city as a planning concept for sustainable urban development is likely to bring up Flanders as a counterexample of what compact urban development represents. On a European ranking of most fragmented countries, Belgium is listed second, after Luxembourg (Jaeger, 2011). To understand the historical development of urban sprawl in Flanders, its current state, and its persistence, De Decker (2011) provides a comprehensive review of spatial planning measures imposed over the last decades -or the lack of them, and their impact on urbanisation. A long-term history of fostering commuting and home ownership and the promotion of a so called 'antiurban' model, aimed at keeping the working class outside the city, seem to have rooted urban sprawl as Flemish cultural heritage. Within Flanders strong spatio-temporal differences in built-up density and morphology can be observed. On a regional scale, high density ribbon development characterizes the central part of Flanders, whereas scattered development of built-up area typically occurs in the province of West-Flanders and in the Campine region east of Antwerp (Verbeek et al., 2010). A recent study of Tempels et al. (2011) examines Flanders' morphological dynamics of urbanization for a selection of rural municipalities, considering built-up density and built form (ribbon vs. scattered development). With regard to built-up density, the position within the Flemish urban network (suburban vs. peripheral areas) and the applicable building regulations during the time period in which the development took place are found to be the main drivers of densification.
2 Driven by a growing awareness of the importance of monitoring urban development, a number of studies have quantified current urban extent and its past growth. However, different definitions of built-up area, different data sets used, other reference years and spatial resolutions hamper a proper comparison of the results obtained in these studies and lead to a certain confusion as well as to some misconceptions in the public debate on urban growth in Flanders and its physical manifestation. In most studies urbanised space refers to areas in which built-up land use is dominantly present (mainly residential, industrial and service-related activities) and which, as such, also includes unbuilt areas that are part of this function-based zoning (like e.g. green spaces in residential areas). This way of looking at urban development contrasts with the mapping and monitoring of areas that are covered by built-up parcels, or areas where actual soil sealing occurs as a result of the presence of dwellings, transport infrastructure, parking lots, industrial and commercial buildings, etc. Depending on the approach taken very different total coverage and growth figures may be obtained that are not directly comparable. Which approach is most relevant clearly depends on the objectives of the study.

3
De Meyer et al. (2011) analyzed the evolution of gross built-up area, defined as the number of 1ha grid cells containing one or more building(s), based on a 12% subset of 1ha grid cells covering the entire Flemish region, and using ortho-imagery as data source. As the study focused on the evolution of built-up area only, transportation infrastructure was not considered. Following this definition, more than one third (36.9%) of the studied area was labelled as built-up, and an average increase of 3.20% since 1990 was observed. Whether or not this number can be considered as representative for the entire region is questionable, as no major urban centres are covered by the sample. Further on, this number solely represents the gross areal expansion of built-up areas, while densification of grid cells which previously contained at least one building is not taken into account. According to the FPS Economy (FOD Economie, 2011), built-up parcels occupied 26.6% of Flanders and 78.9% of the Brussels Capital Region (BCR) in 2011 (27.1% overall), with a respective growth of 35. 6% and 5.8% observed since 1985 (34.3% overall). Although this data set allows monitoring of trends at the regional scale for a period of nearly 30 years, disaggregated data at the level of municipalities is only available since 2003.

4
Following the definition of artificial surfaces in the CORINE land-cover map of 2006, which includes the urban fabric, industrial, commercial and transport units, mine, dump and construction sites and artificial, non-agricultural vegetated areas, 26.9% of the Flanders region and 86.6% of the BCR is covered by these land-use categories (27.6% overall). Furthermore, a respective growth of 5.8% and 1.9% since 1990 is observed for Flanders and the BCR (5.7% overall). For this estimation, it should be mentioned that a minimum mapping unit size of 25 hectares is imposed, resulting in small urban areas and parts of the characteristic ribbon development not being identified as urban, as is confirmed by visual examination. A spatially and thematically more detailed land-use map (10m resolution; 90 distinct classes at the most detailed level, of which 30 are urban), developed by the Flemish Institute for Technological Research (VITO), identifies 30.9% as urban space (i.e. all non-natural and non-agricultural land uses) in 2014 for Flanders and the BCR (Poelmans & Van Daele, 2014). In this case, however, growth figures are not available due to the lack of data required to produce such detailed land-use map for earlier moments in time. Despite providing insight on the spatial extent and expansion of the built-up area or of urban activity related land uses, none of the above mentioned references offers information on the intensity of land consumption in terms of soil sealing.

5
With the emergence of soil sealing fraction as an important environmental indicator in the mid-90s (Arnold & Gibbons, 1996), and the increased availability of high-resolution remote sensing data, there is a growing interest in the monitoring of soil sealing worldwide (Weng, 2012). Besides the consumption of open areas, a less efficient use of natural resources, and fragmentation and reduction of habitat space, soil sealing is detrimental to the soil's environmental functions, e.g. safeguarding soil biodiversity and functioning as a carbon sink (EEA, 2006;Scalenghe & Marsan, 2009). Soil sealing also has a direct impact on the urban climate, as a strong linear relationship between soil sealing and the surface urban heat island can be observed (Yuan & Bauer, 2007). Soil sealing affects the water balance by generating increased runoff (Ungaro et al., 2014) and consequently increased non-point source pollution, and by obstruction of water infiltration and recharge of groundwater bodies, which contributes to groundwater drought (Dougherty et al., 2004;Van Lanen & Peters, 2000). To assess the impact of urbanization at the local scale, detailed information on the density and distribution of soil sealing is thus required (Dams et al., 2013;Poelmans et al., 2010). 6 In order to gain insight in the total amount of sealed surfaces and their spatial distribution, De Meyer et al. (2011) produced a soil sealing map for Flanders and the BCR. Topographic maps (1:10000) were rescaled to a binary (sealed vs. natural) raster (1 meter resolution) and subsequently aggregated to cells of 1 km. Based on field measurements and manual digitizing of sealed areas using ortho-imagery, underestimations of sealed surface cover, caused by small structures such as driveways, terraces, garden houses, etc. not being represented on topographic maps, were accounted for by applying an overall correction factor. Following this approach, 12.9% of Flanders and 63.7% of the BCR was estimated as being sealed (13.5% overall). However, these results have to be interpreted with some caution as (1) the correction factor was determined using a sample of the studied area and applied overall, and (2) due to the time range in which the consulted topographic maps were produced (not all maps refer to the same moment in time) new built-up structures may have been missed, which consequently also affects the correction factor that has been applied.

7
As an alternative for topographic maps, high-resolution Earth observation data, with spatial detail ranging from a few meters to sub-meter resolution, can be used to produce sealed surface maps. Generally, two image classification approaches can be applied. Pixelbased classification mainly relies on spectral information to assign a unique land-cover type to each pixel. To deal with the typically speckled appearance of per-pixel image classifications of urban areas, which is caused by the strong spectral/spatial heterogeneity of urban scenes, rule-based post-classification techniques can be applied, removing the speckle and improving overall classification accuracy, as demonstrated for (a part of) the city of Ghent by Van de Voorde et al. (2007). Alternatively, object-based classification approaches can be applied. Here the image is first segmented into spectrally homogeneous objects, after which each object is assigned to a single land-cover class, taking spectral as well as other criteria, like texture, size and shape of the object, into account. While the object-based approach has been reported as effective for the mapping of scenes that are characterised by a strong spatial and spectral heterogeneity, the major bottleneck of object-based classification is the quality of the segmentation, for which unfortunately no standard approach is available (Neubert et al., 2008). Comparing pixelbased and object-based classification for an urban area, using Quickbird imagery, Yuan & Bauer (2006) obtained a slightly higher overall accuracy and a spatially more homogeneous output using an object-based classification approach. On the other hand, Thomas et al. (2003) demonstrated that applying a rule-based post-classification approach on the outcome of a high-resolution, pixel-based classification, may actually lead to a better classification result than applying an object-based classification approach.

8
Despite the potential of high-resolution image data to map sealed surfaces at a high level of spatial detail, the limited footprint, high cost and time intensive processing of such images hampers their use at a regional or nationwide scale. In addition, the limited historical archive of high-resolution image data restricts their use for spatio-temporal monitoring. For mapping sealed surface cover for larger areas or for studying changes in sealed surface cover over a significant period of time, medium resolution remote sensing data seem therefore more suited. However, in strongly fragmented landscapes, as is the case in an urban and peri-urban environment, the larger pixel size of medium resolution imagery will result in the omnipresence of mixed pixels. The spectral response of such pixels is a combination of the spectral responses of each distinct land-cover type found within the pixel. To deal with these mixed pixels, subpixel classification techniques can be applied, enabling estimation of the fraction of sealed surface cover present within each pixel.
10 Perhaps the most straightforward method to estimate the sealed surface cover fraction from a pixel's spectral properties is the use of linear or non-linear regression (Sawaya et al., 2003;Yang & Liu, 2005;Bauer et al., 2008;Van de Voorde et al., 2011), where the fraction of sealed surface cover or its complementthe vegetation fraction -is directly inferred from the pixel's reflectance in one or more spectral bands, and/or from spectral indices that can be related to the sealed surface or vegetation fraction. Yang & Liu (2005) propose the use of tasselled cap brightness and greenness to estimate the fraction of sealed surface cover from Landsat TM/ETM+ imagery for Pensacola, Florida (US) for two different moments in time to identify hot spots of urban growth. Bauer et al. (2008)  Study area, data and pre-processing of data 14 The Flanders region is situated in the northern part of Belgium, enclosing the centrally located BCR. The region is characterized by a scattered pattern of small, medium and larger cities, well connected by a dense road and railway network, fragmenting the open space. According to the 100m resolution land-use map produced by VITO in the framework of the BELSPO-GroWaDRISK (Drought-related vulnerability and risk assessment of groundwater resources in Belgium) project, nearly one-third of Flanders and the BCR is covered by land uses that are linked to urban activities, with a dominance of residential land use (74%), followed by infrastructure (13%), commercial (9%), industrial (2%) and port areas (2%). Although the Flanders and Brussels regions occupy only 45% of Belgium's total area, with approximately 7.57 million inhabitants in 2014, both regions accommodate about two-third of the total population of Belgium. Flanders has a population density of approximately 472 persons per km², which is higher than the Belgian average (212 persons per km²). Taking into account the large areal extent of the Sonian forest in the BCR, its population density of 7161 persons per km² can be considered very high (FOD Economie, 2014). A high concentration of residents is found in the Flemish diamond, shaped by the central Antwerp-Brussels axis, and expanding to Leuven and Ghent (figure 1). High population densities are also observed north of Hasselt and Genk in the East, around the cities of Kortrijk and Bruges in the West and along the coastline. Since 1981 the total population within the two regions increased by nearly one million (NIS, 1982;FOD Economie, 2014). High population growth is observed in the urban fringe zones, the northern part of the Campine region and along the coastline.
Use of Earth observation for monitoring soil sealing trends in Flanders and B...
Belgeo, 2 | 2016 Figure 1. Study area, population density and location of the IKONOS imagery used as reference data. 15 To obtain a multi-temporal medium resolution coverage of the Flanders and Brussels regions, use was made of the Landsat image archive. Landsat scenes were selected to cover a maximum time span with interval periods of approximately ten years, i.e. 1976, 1987, 2001 and 2013. Prior to the processing, georegistration to the Lambert 72 coordinate system was applied to all image scenes, using a set of ground control points. The characteristics of the selected Landsat scenes, as well as information on the error involved in georeferencing, are documented in table 2.  16 Due to spectral similarities of sealed surfaces and bare soil, for each year in the time series an urban mask was created to avoid potential confusion with bare soil in the estimation of SSF's from the Landsat data. To obtain a base urban mask for the 2013 imagery, an extensive set of ancillary data available for the study area was used to identify those pixels (potentially) containing sealed surfaces. For Flanders, use was made of all artificial feature classes in the large-scale reference data set (GRB-Vlaanderen) which is available for the entire region since 2013. For the BCR, artificial structures were extracted from the UrbIS 2012 large-scale reference data set. Remaining, potentially sealed areas, such as parking lots, driveways, etc., which are not included in the largescale reference data sets, were extracted from a high-resolution (10 meters) land-use map provided by VITO (Van Esch et al., 2011). Finally, an NDVI threshold, separately defined for each Landsat scene, was applied to exclude pixels composed of solely natural surfaces.
17 In order to obtain a consistent sequence of masks for 1976, 1987 and 2001, preceding masks were assumed to be composed of a subset of the more recent mask, i.e. only conversions from natural to built-up area were assumed to have taken place. A reverse two-step filter was applied to identify pixels potentially containing sealed surface cover for previous time steps. Firstly an NDVI threshold, uniquely defined for each Landsat scene, was applied to subtract pixels fully covered with vegetation from the 2013 base mask. Secondly, all pixels more than two standard deviations away from the linear trend between the NDVI values of two overlapping scenes for two consecutive periods in the time sequence were removed. Based on visual examination the latter criterion proved to be able to exclude former agricultural areas from the mask. For the 1976 urban mask, which was produced at a resolution of 60 meters in order to match the spatial resolution of the Landsat 1 and Landsat 2 imagery, NDVI values for pixels of the 30m resolution urban mask of 1987 were averaged to define the linear NDVI relationship with the corresponding 60m pixels in the 1976 image. The four urban masks that were obtained, in which a complete absence of bare soil is assumed, thus comprise all pixels (partly) containing sealed surfaces, and delineate the area for which SSF's for each year in the time series are estimated. binary maps were then spatially aggregated to match the spatial resolution of the Landsat imagery, whereby average vegetation fractions and NDVI-values were assigned to each pixel. This procedure is illustrated for a subset of the Ghent reference area in figure 2. Figure 2. Extraction of reference vegetation fractions from Ikonos imagery using ancillary data and NDVI-thresholding. 19 A stratified random sampling was applied to extract an independent training and validation set from the Landsat pixels being part of the urban mask for which reference vegetation fractions were known. Given the difference in acquisition dates of the Ikonos images (acquired in the 2002/2003 growing season) and the Landsat data, and consequently the potential discrepancies in vegetation cover between both image sets due to land-cover changes or seasonal differences, a selection criterion was defined to filter out pixels from the training and validation set for which the spectral properties had changed between the acquisition dates of the Landsat and the Ikonos imagery. More in particular, all pixels more than two standard deviations away from the linear trend between the Ikonos NDVI and Landsat NDVI values were removed from the sample ( Van de Voorde et al., 2008). This procedure results in an independent set of training and validation data for unmixing each Landsat scene.
Linear regression unmixing 20 For each Landsat scene, a linear regression analysis was performed using the known vegetation fraction derived from the IKONOS imagery as dependent variable and the Landsat NDVI as independent variable: 21 with VF i the reference vegetation fraction in pixel i, a 0 the model's intercept, a 1 the regression coefficient, NDVI i the NDVI-value of pixel i and ɛ i the residual for pixel i. The inclusion of spectral bands in the regression equation, and, as such, the extension of the simple linear regression model into a multiple linear regression model, resulted in little to no significant improvement in the estimation of the vegetation fraction, depending on the Landsat scene considered. Hence, for reasons of simplicity and consistency, for each Landsat scene a linear regression model including the NDVI as the sole variable was defined and consequently applied to unmix each image. 22 To assess the performance of the scene-dependent regression models on the independent validation data, the mean fractional error (ME) and the mean absolute fractional error (MAE) were calculated.
Use of Earth observation for monitoring soil sealing trends in Flanders and B... Belgeo, 2 | 2016 23 with N the total number of validation samples, VF i the estimated vegetation fraction for pixel i, and VF' i the observed vegetation fraction for pixel i. 24 The ME provides a measure of bias in the fraction estimation. Hence, the obtained residuals are expected to be normally distributed, with an ME close to 0, exhibiting no systematic under-or overestimation. The magnitude of the estimation error is assessed by the MAE. In order to provide insight on the error distribution of the fractions for each image, the standard deviation with respect to the ME is documented as well. Since ME values prove to be very small (see below), indicating that almost no bias is present in the results, standard deviation will be very similar to an RMSE estimate of the error, attributing more weight to larger errors than the MAE. Finally, the correspondence between observed and predicted vegetation fractions is expressed by the Pearson correlation coefficient (R).
25 Following the earlier stated assumption that each pixel within the urban mask is composed of vegetated and/or sealed surfaces only, the estimated SSF was obtained by taking the complement of the vegetation fraction: Temporal consistency 26 To produce a consistent sequence of SSF maps, a levelling was applied to those pixels in which the estimated SSF in 2001 and 2013 was lower than the fraction estimated in the preceding maps of 1987 and 2001, respectively. Analogue to the delineation of the historic urban masks, it was assumed that the transformation of natural to sealed surfaces is the only possible land-cover change. Although occasionally the opposite transformation might occur, it is relatively rare. As such a decline in SSF at pixel level from one period to the next will be more likely to be caused by error in SSF estimation at the level of the individual pixel than being the result of an actual decline in SSF. The assumption of an increase in sealed surface cover as the sole direction of change was considered as the most appropriate strategy for analyzing spatio-temporal changes in sealed surface cover at the regional level. Given the lower spatial resolution of the 1976 images, the choice was made to apply a reversed levelling of SSF estimates for this date, i.e. the aggregated SSF estimations of 1987 were used to define the upper limit for the estimated SSF in 1976.

Results and discussion
Urban mask delineation 27 Using ancillary data to delineate the current state urban mask results in a masked area of 5109 km² (37.13%), potentially containing sealed surfaces. Applying the reversed two-step filter approach to obtain urban masks for previous time steps, total areas of 4407 km² and 4923 km² are obtained for 1987 and 2001 respectively (table 3). This corresponds to an The higher figure obtained for the 1976 urban mask can be attributed to the lower spatial resolution of the Landsat 1 and Landsat 2 imagery, which obviously results in a larger area covered by cells potentially containing sealed surfaces. Table 3. Sealed surface area, density, growth and share, globally and for different land uses.
28 Because the unmixing model used in this study assumes urban mask pixels to be composed of vegetated and/or sealed surfaces only, the extraction of a representative urban mask for each time step is of utmost importance. Figure 3 demonstrates the application of the proposed approach for delineating urban masks for each time step, starting from the 2013 urban mask. For 2001 and 1987, the shrinkage of the urban mask while going back in time reflects well the urban expansion and the filling-up of existing urban areas, with linear structures present over the entire period clearly being shown at all time steps. Due to the coarser resolution of the Landsat MSS imagery, the 1976 urban mask shows clear differences in urban pattern compared to the other three masks: large urban areas and built-up structures appear 'blown-up', whereas small linear objects are more scattered. A quantitative validation of the urban masks for different time steps is difficult to achieve, given the absence of detailed historic data at the regional level. Visual validation of the urban masks for selected areas, however, demonstrates the applicability of the approach proposed. Nevertheless, potential errors in the delineation of the maskswhich will inevitably be presentask for caution when using the results obtained in this study for analysing changes in sealed surface cover through time at the local level. This is particularly the case for changes reported between 1976 and 1987.  30 Several factors may explain the observed error in SSF estimation, yet the most important factor in regression-based approaches contributing to uncertainty in the modelling is likely to be the geometric mismatch between the Landsat image data and the dataset from which the reference SSF fractions are derived, in this case high-resolution Ikonos imagery. As is illustrated by the georegistration errors obtained for each Landsat image (table 2), such mismatches are impossible to avoid. Error in the georeferencing of the imagery will add noise to the regression model, resulting in more dispersion along the trend line, yet it will also introduce uncertainty in the validation of the obtained models at pixel level, where errors observed will be partly caused by model error, yet will also partly result from the geometric mismatch between image and reference data. In other words, error reported at pixel level is likely to be inflated by this effect. 31 In order to obtain a better idea of the uncertainty involved in SSF modelling and validation at pixel level, for each pairwise combination of overlapping Landsat images available for the same time step, a random sample of 1000 pixels, part of the urban mask, was selected. The distribution of absolute differences in sealed surface fraction estimates for overlapping pixels is shown in figure 4, for each time step. As can be seen, the distributions show median values around .20 and a distance between the first and third interquartile of about .30, with a stronger concentration on the lower side of this range. These differences in reported error at pixel level are to a large extent caused by error in the georeferencing of the Landsat imagery and will, in our study, also introduce uncertainty in the temporal changes observed at pixel level.  Figure 4. Distribution of absolute differences in sealed surface fraction estimates for pixels in the overlapping area of two Landsat images, for each time step. 32 The problem of the geometric mismatch between image and reference data has been reported in many unmixing studies, and explains why validation of unmixing is often done for larger cell sizes (Song, 2004), or at the level of meaningful spatial units such as urban blocks or neighbourhoods (Van de Voorde et al., 2008;Okujeni et al., 2014), by aggregating the fraction values obtained at pixel level. This way the impact of coregistration error on the validation is reduced. However, this still does not solve the problem of the geometric mismatch in the calibration phase of regression-based modelling. One way to deal with this mismatch is to calibrate regression-based unmixing models based on synthetic mixture data, as was recently suggested by Okujeni et al. (2014). As this approach no longer requires the use of a reference fraction dataset in the calibration phase of the modelling, it also effectively deals with the problem of the geometric mismatch.
33 Apart from the impact of geometric shifts, other factors that introduce uncertainty in the modelling approach and in the trends observed at pixel level may relate to sensor characteristics (point spread function, IFOV geometry), seasonal effects (phenology of the vegetation), accuracy of the reference fractions (derived from the Ikonos imagery), and inconsistencies in the urban mask definition over time. Although the influence of each of these factors on the accuracy of the unmixing is an interesting research question on its own right, it is beyond the scope of this paper. In a study on the estimation of the vegetation fraction from Landsat ETM+ data for the Brussels Capital Region, using linear regression unmixing in a similar way as in the present paper, Van de Voorde et al. (2008) demonstrated that aggregating the obtained fractions from the pixel level up to the level of urban neighbourhoods reduced the mean absolute fractional error from .10 to .02. As in this study the evolution of soil sealing is analysed at an aggregated level (regional,  36 To validate the entire processing chain proposed in this study in a more rigorous way, it would be interesting to confront the obtained SSF maps with independent data on the presence of sealed surface cover not obtained through remote sensing. However, as none of the large-scale reference datasets that are currently available take actual soil sealing into account, the only way to perform a validation independent of the procedure used in this study would be through manual interpretation of high-resolution aerial photography for a sufficiently large set of sample areas representing different built-up densities and patterns of settlement for each step in the time series. Apart from the difficulty of finding appropriate data to perform such a validation, it would constitute a huge effort, for which the resources that would be required are currently not available. Land-use specific sealed surface trends 37 When considering the relative share of each land-use class in the total area covered by sealed surfaces, residential areas constitute the largest fraction (47.5% in 2013), although their share has been slightly decreasing since 1976 (table 3). The contribution of infrastructure, commercial, industrial and port areas to sealed surface cover within the two regions has remained more or less stable (21%). Sealed surfaces also occur within other land-use classes (agriculture, forest, parks, etc.). These sealed areas mostly correspond to local streets and their adjacent buildings (82% of these sealed surfaces are located within a distance of 100 meters from the street network), which due to their low built-up density, in combination with the rules applied to produce the land-use map, are located outside the five higher mentioned built-up land uses. Within the framework of Use of Earth observation for monitoring soil sealing trends in Flanders and B... Belgeo, 2 | 2016 this paper, all sealed surface area outside the built-up land-use classes was therefore considered to be part of a mixed land-use class. 38 Analyzing the presence and growth of sealed surfaces for different land uses shows some clear deviations from the overall sealed surface density and growth trend. Between 1976 and 2013, the largest absolute increase of sealed surfaces occurs within residential (522.5 km²), mixed (447.7 km²) and commercial (109.1 km²) areas, whereas the strongest relative increase is observed for mixed (+112.1%), commercial (+99.3%) and industrial (+75.1%) areas. Within residential, infrastructure and mixed land uses a decreasing growth rate is noted. Within employment related land-use classes (commercial, infrastructure and port related activities) an initial increase in growth rate is observed between 1976 and 2001, after which a decrease sets in. These trends should however be interpreted with caution, as the land-use map used in our analysis solely describes the current spatial extent of each land-use class. Consequently, a distinction between sealed surface growth linked to a functional change from one built-up land use to another or to densification of a particular land use in areas where it was already present before cannot be made. A clear spatial pattern is observed when examining the share of sealed surface growth caused by industrial, commercial and port related activities at the level of the municipalities. The highest shares are observed in Antwerp and Bruges (port expansion), as well as northeast of Brussels and Kortrijk (commercial development). High shares are also found in the Flemish diamond, around Ghent (port expansion) and along the Albert canal (figure 7c).
39 Although roughly one fifth of the total study area is sealed, the sealed surface density within built-up land-use classes is much higher. Highest average densities are obtained for port areas (81.6%), industrial (70.2%) and commercial (59.9%) areas (table 3). The probability density functions of SSF's within each land-use class show a nearly normal distribution of SSF within residential and infrastructure areas, whereas within port and mixed areas a dominance of respectively high and low fractions appears (figure 8). For commercial and industrial areas a steady increase in the frequency is noted for increasing SSF values. These contrasting patterns point out inter-class (e.g. large built-up units in industrial areas compared to smaller detached residential dwellings) as well as intra-class differences in SFF (e.g. low versus high dwelling densities) which can be linked to typical built-up morphologies. Keeping in mind earlier remarks on the lack of historic land-use data, historic sealed surface densities for different land uses reported in table 3 should be interpreted with caution. It would nevertheless be interesting to include detailed historic land-use data in the analysis to better examine whether or not and to which degree soil sealing density has changed within distinct types of land use.  1976,1987,2001 and 2013 have been used in combination with high-resolution Ikonos data to estimate subpixel sealed surface fractions for different time steps. Based on the validation of the unmixing models and visual examination of the obtained sealed surface fraction maps, the linear regression unmixing approach proposed can be considered a simple and reliable method for mapping and monitoring of sealed surface fractions at a regional scale. It should be noted though that the quality of the urban mask that is used to delineate areas on which the unmixing is applied is of utmost importance in the outlined strategy, and largely depends on the availability of ancillary data identifying potentially sealed areas for the last step in the time series.
41 For 2013 an overall sealed surface cover of 19.5% of the area of Flanders and the Brussels Capital Region is obtained, corresponding to a total sealed surface area of approximately 2687 km². Since 1976, the sealed area has increased by 82%, representing an extra 1214 km² of sealed surface cover. Nevertheless, a decreasing growth rate, which currently amounts to 6.10 hectares per day, is observed. Residential areas account for nearly half of all sealed surface cover, and have shown the largest absolute increase over the past 37 years. Sealed surface density varies greatly depending on land use and is the highest for port related land uses (81.6%), followed by industrial areas (70.2%). Within Flanders, a large spatial variability in sealed surface density, as well as in the pattern of sealed surface growth is observed. The strongest growth in sealed surface cover is observed in the fringe zones of urban agglomerations and in strongly industrialized regions.
42 It would be interesting to analyse the sealed surface fraction maps obtained in combination with detailed data on the characteristics of built-up areas (e.g. year of construction and built-up morphology), population density and spatial context. This might offer more insight in the relationship between sealed surface density and different urban typologies, and in how changes in sealed surface density and in growth patterns through time can be explained. Such insights could be used to estimate future sealed surface growth for different planning scenarios and to assess related environmental impacts, e.g. impact on surface water runoff and on recharge of groundwater bodies.

ABSTRACTS
The on-going growth of urban area in Flanders and in the Brussels Capital Region over the past decades has resulted in a highly sprawled urban tissue, consisting of large and smaller urban agglomerations, connected by a well-developed transportation network. The conversion of open land to urban area is accompanied by an increase in soil sealing, affecting the hydrological cycle and the urban climate. Despite a growing interest in monitoring the process of soil sealing in urban areas, to date no detailed information on the presence and evolution of sealed surfaces is available for Flanders. In this paper a linear regression unmixing approach is proposed to map and monitor changes of sealed surface cover at the regional scale, using medium as well as high resolution remote sensing data. Applied to Flanders and the Brussels Capital Region, a total sealed area of 2687 km² for 2013 is found, corresponding to an increase of 82% since 1976. Residential areas account for nearly half of the sealed area and show the largest increase in sealed surface cover over the past 37 years.