Geoffrey Ginoux
email: geoffrey.ginoux@univ-reims.fr

Isabelle Vroman

Sébastien Alix

INFLUENCE OF FUSED FILAMENT FABRICATION PARAMETERS ON TENSILE PROPERTIES OF POLYLACTIDE/LAYERED SILICATE NANOCOMPOSITE USING RESPONSE SURFACE METHODOLOGY

This study aimed at assessing and optimising the influence of printing speed and extrusion temperature in a fused filament fabrication (FFF) process on the tensile properties of a polylactide/layered silicate nanocomposite. Mathematical models using Doehlert designs were formulated to examine factor and interaction effects. The models were corroborated by measurements using capillary rheology, tomographic images, and crystallinity analyses to find physical explanations for the differences in tensile properties. The tensile properties were a non-monotonic function of printing speed, which may be due to various deposition defects that influence the porosity of composite tensile specimens. This study provides new insights into FFF process optimisation regarding rheological behaviour and mesostructure of nanocomposite by highlighting new modes of deposition defects that originate from process parameter settings and materials. The results contribute to the properties mastery of FFF-processed materials.

Introduction

Additive manufacturing (AM) techniques are currently the subject of considerable interest due to their promising applications. Initially used for rapid prototyping, AM makes it possible to manufacture complex parts because its layerby-layer approach enables the manufacture of expensive end-use products at a lower cost or of parts impossible to manufacture using traditional subtractive or formative techniques. While AM refers to several technologies 1,2 , fused filament fabrication (FFF) is one of the most commonly used due to its speed and cost-effectiveness 3 . However, additively manufactured parts have inferior mechanical performance compared to those produced by subtractive and formative processes because of manufacturing defects 4 , which hindered the development of this technology for smallscale production or medical applications 5 . Manufacturing defects in those parts include porosity 6 and distortion 7 , which can initiate fracture during mechanical tests. Those defects are due to the spreading of the molten material in successive layers which causes an anisotropic structuring, provoked by the raster angle during the process 8 . Currently, an important topic in this field is predicting the mechanical properties as a function of printing parameters to limit those defects and to obtain desirable properties. Numerous parameters can be adjusted to modify the mechanical properties of the product, including shell number, layer thickness, infill percentage, deposition angle, build orientation, raster width, displacement speed of printing head, and temperature of liquefier 9 . Those parameters modify the structure of the built parts, namely porosity and crystallinity 10 , and thus the mechanical properties 11 . Printing parameters may influence the thermal cycles and cause a differential of shrinkage, which could cause distortion, internal residual stress and lower mechanical properties 12 . Printing speed and extrusion temperature have been shown to impact the adhesion between deposited filaments 13 , and are suspected to indirectly modify the structure by influencing the ability of the filament to be correctly deposited 14,15 . To date, the impact of printing speed on mechanical properties is not clear, and has been shown to have both positive and negative effects [14][15][16][17][18][19] . In general, it seems that a low printing speed leads to better adhesion between layers. This effect was reported by Deng et al. 14 . They reported that low speeds lead to uneven deposition. Other studies showed a better adhesion at low speeds 18,19 . Sukindar and Jaafar 15 stated that high speeds reduce the consistency of the extrusion process, which would lead to turbulent flow and irregular deposition. Regarding extrusion temperature, Spoerk et al. 20 reported a positive effect, since an increase in temperature decreases viscosity and promotes filament adhesion.

The poor filament adhesion would increase the anisotropy of the printed part by lowering the mechanical properties in the building direction 21 . Other phenomena might also influence the effect of temperature on mechanical properties, such as crystallinity 10,22 or matrix degradation 23 , and cause a more difficult understanding of this thermal behaviour. Printing speed and temperature are thought to interact, since printing speed can modify the time between two adjacent depositions, which is the time allowed for a filament to weld before the temperature falls below the glass transition 24 . The dependence of mechanical properties on thermo-kinetic parameters is suspected to be more nuanced than expected, and there may be an optimum combination of properties.

Design of experiments (DOE) is a technique that is often used to model the effects of FFF process parameters and their interactions both cost-and time-efficiently. The response surface methodology (RSM) is particularly effective for observing an optimum response. The most well-known DOEs in the field of additive manufacturing are central composite design, the Taguchi method, and the Box-Behnken method 25 . Doehlert design 26 is also interesting due to its uniform space-filling and iterative nature, which makes this design effective to find optima, especially in case of complex effects weakly predictable with a single quadratic model. This iterative fashion allows to start with a restricted domain and to widen it only when necessary, improving the time and cost efficiencies of the experiments when the localisation of optima is multiple, i.e. visible in several domains, or not obvious. Surprisingly, while this method is not commonly used in AM studies, it is well known in other fields like chromatography 27 . Existing models are essentially focused on the mechanical optimisation of neat polymers, while nano-filled polymers can also be used to compensate for the previously cited loss of mechanical performance 2,28 , and to help functionalise the materials 29 or ease the processability 30 . However, filled polymers are known to be stiffer than neat polymer, especially while additively manufactured, and need to compensate the lack of ductility with an optimisation of the process 31 . A Doehlert design was then chosen to understand how printing speed and extrusion temperature influence the mechanical properties of printed polylactide (PLA) filled with layered silicates. Rheology is dependent on those factors, which could modify the structure and thus the tensile properties. Hence, capillary rheology, X-ray micro-computed tomography (microCT) and differential scanning calorimetry (DSC) analyses were used to propose physical explanations for the mathematical models obtained with DOE.

Materials and methods

Materials and nanocomposite preparation

PLA was purchased from NaturePlast (France) in pellet form as extrusion-grade PLE 005. Table 1 summarises some its properties. PLA was first dried at 70 °C overnight before nanocomposite compounding with a co-rotating twin-screw extruder (BC21, Clextral, France; L = 900 mm, L/D = 36) with temperatures set between 165 °C and 180 °C at 300 rpm. The residence time was approximately three minutes with a mass flow of 1.8 kg.h -1 . A masterbatch of clay at 13.2 wt% in PLA was obtained. The masterbatch was then diluted with neat PLA in the same extrusion process at a desired loading of around 3 wt%, which is an optimal loading for thermal and mechanical properties 32,33 . Nanocomposite was then pelletised, and the MMT mass fraction was controlled by thermogravimetric analysis for residue after degradation at 800 °C for 10 min under air flow. MMT mass fraction was determined to be 2.7 ± 0.2 wt%.

Filament fabrication

Nanocomposite pellets were first dried overnight at 70 °C in a ventilated oven and then shaped into a 1.65 ± 0.05 mmdiameter filament using a single-screw extruder (BX10, Axon Plastics Machinery ab, Nyvang, Sweden; L/D = 18) with a screw speed at 20 % of maximum and a 460 g/h mass flow. The profile temperature, separated into three zones, was set to 190 °C, 185 °C, and 185 °C from the feeding zone to the nozzle. The filament was pulled by a Haake™ PP1 (ThermoFisher Scientific, Waltham, MA) to obtain the required diameter.

Fabrication of specimens

Once the filament was obtained, tensile test specimens were manufactured using a Replicator 2 (MakerBot Industries, LLC, Brooklyn, NY) with a 0.4 mm-diameter nozzle. The specimens were first modelled using Solidworks ® software (Dassault Systèmes, SolidWorks Corporation, Waltham, MA) with dimensions according to 5A-type specimens from ISO 527-2 standard. The STL file was then imported to Makerbot Print software. Specimens were printed with two shells, in a 90°/0° crossed-deposition fashion (see Figure 1), which gives greater modulus of elasticity than 45°/-45° deposition 31,34 , with a 0.2 mm layer-height, for a total of ten layers. Printing speed has been modified according to the experimental design but is constant during printing, i.e. it is identical for all the layers of the sample. Build plate temperature was at room temperature. Infill percentage was set to a theoretical 100 %.

Δ𝐻 𝑚 is the measured melting enthalpy, Δ𝐻 𝑐𝑐 is the measured cold crystallisation enthalpy, Δ𝐻 𝑚 0 is the melting enthalpy of a theoretical 100 % crystalline PLA (91 J.g -1 ) 35 , and 𝛷 𝑚 is the mass fraction of MMT.

Doehlert design

The Doehlert method is an iterative technique for response optimisation. For this study, one iteration was performed. The case with two factors (k) uses a simplex triangle and requires seven experiments (N = k² + k +1). Since the method is iterative, the values must essentially consider the "size" of the simplex (geometric unit). Differences of 15 mm.s -1 for printing speed and 25 °C for extrusion temperature were a good compromise between a sufficiently large gap to observe a difference and a not-too-large gap to avoid missing a local optimum. The iteration included four experiments, from the first design to three new ones (Table 2 and Figure 2).

Table 2. Process parameter values with test numbers for both Doehlert designs and respective configuration numbers for each design, with printing speed as factor 1 and extrusion temperature as factor 2. All calculations were performed with Excel ® software (Microsoft Corp., Redmond, WA) using matrix calculations. For each design, the model was written as

Test

𝑦 𝑡ℎ = 𝑎 0 + ∑ 𝑎 𝑖 𝑥 𝑖 𝑘 𝑖=1 + ∑ 𝑎 𝑖𝑖 𝑥 𝑖 2 𝑘 𝑖=1 + ∑ ∑ 𝑎 𝑖𝑗 𝑥 𝑖 𝑥 𝑗 𝑘 𝑗=𝑖+1 𝑘-1 𝑖=1 + 𝜀 (2)
yth is the theoretical response of the experiment at the chosen levels, a0 is the response of the central point, ai and aii are the coefficients of the mathematical model from factor i, aij are the coefficients of the interaction between factors i and j, xi and xj are the factor levels, k is the number of factors, which is 2 in this study, and ε is the residue.

Analysis of variance (ANOVA) was used to assess the significance of the coefficients and models and normality of the results were checked using Shapiro-Wilk test.

Results and discussion

Statistical analysis

Tensile test results for the Doehlert designs are summarised in Table 3. Young's modulus (E) and ultimate tensile strength (σM) were chosen as responses for the models. Elongation at break and elongation at ultimate tensile strength were not treated by the Doehlert method due to high experimental standard deviations, leading to low significance for all Printing speed (coded units)

the coefficients and models. Data for both designs and both responses follow a normal distribution according to Shapiro-Wilk test at a 0.05 level. A1 to B2). The models with ultimate tensile strength as a response for both designs are not significant at the 90 % confidence level (Tables A2 andB2). The low significance could be attributed to the unsuitability of the model for predicting the variation of ultimate tensile strength with printing speed and extrusion temperature, or to the high experimental standard deviations compared to the variation of the response between each configuration. On the contrary, models with Young's modulus responses are significant at more than the 90 % confidence level. A single optimum was expected from among the disparate results in the literature, but several optima were observed instead (Figure 3a and3b). The first mathematical model shows an optimum, which is a minimum (Figure 3a). To compensate for the loss of mechanical properties with FFF compared to formative processes, a maximum was desired. The model considers a significant positive first-order coefficient for printing speed, a1, while that for extrusion temperature, a2, is not significant (Table A1). This maximum can be expected at higher speeds. The second design was constructed by sliding the first design in an iso-thermal fashion (a2 = 0) to higher speeds (a1 > 0). Three new points were added, while keeping four points of the first design, and matrix calculations were carried out again. The results after this iteration led to the presence of another optimum, which was a maximum. Regarding the coefficient significance for p-values lower than 10 %, the effects of printing speed and extrusion temperature on Young's modulus follow a second-order polynomial for each design (Tables A1 andB1). Young's modulus is a function of printing speed with first-and second-order terms, a1 and a11. The dependence with temperature follows only a second-order term, a22, while the interaction, a12, is null. The difference between the two designs is reflected in the signs of the second-order terms, a11 and a22: they are both positive in the first model and negative in the second model. Those values mean that Young's modulus tends to increase globally with printing speed, but not with temperature or printing speed/extrusion temperature interaction. However, the effect of printing speed is locally more nuanced as several local optima are noticed on each design. The change of sign for both second-order terms, a11 and a22, between the two designs is the reason why a minimum is observed on one design and a maximum is observed on the other one. Although the interaction term is null for both models, the effect of extrusion temperature is dependent on the range of printing speed. Those observations point out a more complex model for predicting Young's modulus plotted against printing speed and extrusion temperature. Indeed, Young's modulus seems to follow globally a third-order polynomial model, which could lead to disparate local effects. As an example, Montero et al. 36 reported that extrusion temperature did not affect tensile strength, whereas Spoerk et al. 20 and Alafaghani and Qattawi 37 stated that extrusion temperature had a remarkable positive effect. Because of the interaction with a large range of printing speed, the effect may be different from one study to another. Similarly, the use of only two values per factor could cause a gap between the actual and observed effects; the presence of an optimum could be hidden and only the effect of the first-order term may be visible. Figures 2a and2b show that printing speed has a negative effect on Young's modulus, followed by positive and negative effects. These different effects of printing speed are in agreement with the literature [14][15][16][17] and are likely dependent on the study conditions. The complexity of the printing speed effect on Young's modulus could be attributed to numerous physical phenomena, such as adhesion, distortion, crystallinity, and so on.

Physical explanation

In the first Doehlert design, the second-order term of the printing speed, a11, had the greatest effect on tensile properties and led to a highly curved surface. Competition between inner residual stresses and bond formation may be important.

Zhang and Chou 38 concluded that increasing print speed led to a greater thermal cycle frequency and higher residual stresses, which worsen mechanical properties. The residual stresses may be due to distortion caused by a gradient of crystallinity 12,39 . DSC measurements of the two first printed layers and the two last layers show different crystallinity rates (Figure 4). The bottom of the tensile test specimen shows a cold crystallisation enthalpy of the same order as the melting enthalpy, while the top shows a similar melting enthalpy, but a lower cold crystallisation enthalpy. As a result, the crystallinity rate is calculated to be around 1.4 ± 1.1 % for the first printed layers and around 38.1 ± 3.3 % for the last layers, which is due to the lower temperature of the build plate compared to the temperature of the deposited layers. The difference of temperature of the subtract leads to a differential of cooling rate, slower for the latest layers, that causes higher crystallinity rate with the height position. A gradient of crystallinity is observed, and this can cause distortion and lower mechanical properties as a function of the printing speed 12 , although composite materials are less subject to distortion than neat materials as shrinkage is less pronounced when materials are filled. On the one hand, a faster printing speed was associated with a higher shear rate. The shear rate was calculated as

𝛾̇= 4𝑄 𝜋𝑅 3 3𝑛 + 1 4𝑛 (3)
Q is the flow rate, which is approximately the product of the printing speed, the layer height and the raster width, R is the radius of the outlet nozzle, and n is the flow behaviour index obtained from rheological measurements (Figure 5). The shear rate in the Doehlert designs was estimated to be between 250 s -1 and 5000 s -1 , which is in the range of measurements obtained by capillary rheology (Figure 5). Thus, an increase in printing speed reduced the viscosity, which would facilitate deposition of the polymer by allowing it to flow more suitably through the nozzle onto the previous layer, similarly to extrusion temperature 20 . Figure 6 illustrates this statement, with each axis view represented by a pair of different configurations. When observed on-edge, especially in the x-and y-axes views where deposition is visible, when printing speed was low (configuration N° 4), deposition defects were apparent. Those defects were due to the polymer not being fluid enough under a low shear rate (Figure 5). Because of those defects, Young's modulus and the ultimate tensile strength tend to decrease at lower printing speed, while the mechanical properties increase with a decrease in residual stresses at a lower printing speed, leading to curvature in the response surface. However, the variation of viscosity with the shear rate did not impact bond formation after deposition, when the shear rate was null. On the other hand, the time between adjacent filament depositions is reduced by greatly increasing the printing speed (up to seven times greater from test N° 4 to test N° 9), which provides less time for the polymer to cool down. These higher temperatures led to lower viscosity, as seen in Figure 5 and in the literature 30,40 . This phenomenon is more relevant for small specimens like 5A-type specimens in which the deposition time between two adjacent filaments is between 0.1 and 5 secondsin this Doehlert networkwhile the time required to reach the glass transition is around 1 to 3 seconds 41 .

Depending on the configuration, a filament could be deposited next to the previous deposition either in a "viscous" or "glassy" state, which, respectively, facilitates or restricts the interdiffusion of polymer chains between the deposited filaments 13 . A more developed adhesion could be expected with configuration N° 1 rather than configuration N° 4, as seen on tomographic images. In this way, increasing printing speed increases the shear rate and reduces the time between two filament depositions and then the loss of heat of the adjacent filaments, which tends to give a more fluid material during and right after deposition and favours bond formation between the deposited filament and the previous and adjacent one.

These observations were inconsistent with those from literature 15,18,19 , who reported that high printing speed is a source of deposition defects and that low printing speed enhances the adhesion between deposited filaments. This difference at low printing speed may come from the ability of the materials to suitably flow over the previous deposited layer, which is dependent on the rheological behaviour of the used materials. A critical viscoelastic behaviour may result from one case to another and defects may or may not appear at lower printing speed with other materials. Differences may particularly be noticed between filled or unfilled polymer as filled polymers tend to have greater viscosity at low shear rate than unfilled polymers. This greater viscosity for filled polymers could overpass the critical viscosity from which materials are not fluid enough to be effectively deposited and/or could prevent the adhesion at zero-shear rate by inhibiting polymer chain interpenetration between adjacent filaments because of the hindering filler 40 . Therefore, defects might result from two different mechanisms with a change in printing speed:

(1) a polymer that is not fluid enough at low shear rate decreases the ability of the filament to be effectively deposited over the whole raster width, especially if materials are filled

(2) deposition that is too fast may lead to turbulent flow and inconsistent deposition.

The interaction between printing speed and extrusion temperature could greatly contribute to this difference in observations. Tsouknidas et al. 16 also observed a monotonous and negative effect of printing speed (from 30 to 220 mm.s -1 ) on the maximum stress. In contrast, the Doehlert designs revealed an initial negative effect, then positive and negative effects from 15 to 105 mm.s -1 . This difference could be due to thermo-kinetic interaction and the use of other printing conditions, since Tsouknidas et al. 16 used a constant temperature of 180 °C for PLA, with an infill percentage of 50 %. Due to this low temperature and partial filling, which accentuates heat loss due to convection and reduces the heat transfer between filaments through conduction 42 , interpenetration of the polymer chains would be limited. Thus, only a phenomenon that gives deposition defects would be involved and would only reduce the mechanical properties as a function of printing speed. The curvature on the surface responses is the expression of several phenomena acting on the deposition of the filaments. The capacity of the filament to be effectively deposited would directly impact the mesostructure and therefore the mechanical properties. This new competition between the ability to be deposited on the whole raster width and the ability to be deposited consistently leads to a new optimum of mechanical properties at higher speeds. Indeed, increasing printing speed would lower the viscosity and facilitate deposition of the material, until the deposition is too fast and causes new defects due to turbulent flow that decreases the adhesion between layers. The optimum would then be a maximum, as observed with the second Doehlert design. Figure 7 illustrates this phenomenon by comparing the fracture zone of the printed tensile test specimens at the same extrusion temperature but different printing speeds. Printing at the lowest speed, 15 mm.s -1 , leads to triangular voids (Figure 7a). Those voids are representative of depositions that do not occur over the whole raster width. The filament is in contact with the previous layer but is not fluid enough to be deposited over the entire raster width. In addition, on the edges of the tensile specimen, the deposition fills gaps with very low precision. At a higher speed, 45 mm.s -1 (Figure 7b), defects are still visible but less pronounced.

Among other features, the contours are better defined, and the overall porosity seems lower. At an even faster speed (Figure 7c), deposition is ideal, with little porosity. The deposition is regular, and the printing accuracy is greatly improved compared to the two previous speeds. This configuration is test N° 1, 75 mm.s -1 and 215 °C, which experimentally gave the highest values. The highest speed, 105 mm.s -1 (Figure 7d), led to the creation of new defects.

The pores outside the contours are not triangular, but rather are elongated, which is attributed to the lack of adhesion of the filament under rapid deposition. This phenomenon does not allow effective contact between the inter-layer filaments.

At high speed, the filament would not interfere with the previous layer over a portion of the raster length. In this way, the porosity increases again, and this explains the inferior mechanical properties of the configuration at 105 mm.s -1 compared to the configuration at 75 mm.s -1 . The two deposition defect modes act on different steps of adhesion 13 . Sun et al. 24 essentially focused on diffusion at the interface in terms of temperature, but modification of printing speed would also act on the bonding quality between filaments via surface contact and the initial neck growth. On one hand, at low print speed, the deposition is in contact, but only partially, with the raster width, which decreases the ability of the deposition to fully develop contact between filaments by reducing the initial neck growth. On the other hand, at high speed, the deposition is inconsistent, with no contact on some portions of the raster, which prevents the adhesion process. The ability of the filament to be effectively deposited would directly modify the porosity of the product. To verify the impact of porosity on mechanical properties, the filling rate was measured. All configurations were tested, and the filling rate was used with the two Doehlert designs. Young's modulus and the ultimate tensile strength are weakly correlated to the filling rate with coefficients of determination equal to 0.3312 and 0.0791, respectively, via linear regression. By applying the RSM with the filling rate in response, it appears that the surfaces show a minimax for the first design and a maximum for the second design at a more than 90 % confidence level (Figure 8). The presence of a maximum could indicate that porosity significantly impacts the mechanical properties, but its effect is counteracted by other phenomena that reduce the correlation between porosity and mechanical properties. Those phenomena can contribute as much to the creation of internal residual stresses as the development of adhesion between deposited filaments. This last phenomenon not only considers the contact surface, depending on the porosity, but also the interpenetration of the polymer chains between the deposited filaments, which is not visible in microCT images. The minimum mechanical properties could be due to competition between internal stresses and fluidification of the material, which are properties at the molecular scale and at the scale of the mesostructure. Therefore, the filling rate may not necessarily pass through a minimum. The porosity could play a role in the optimisation of mechanical properties, but its effect is not dominant over the effects of other phenomena.

(a) (b) The crystallinity rate was measured for all configurations to assess its effect on tensile properties compared to the effect of the filling rate. Wittbrodt and Pearce 43 previously observed a lack of correlation between the crystallinity rate and yield stress with changes in extrusion temperature, but the correlation could differ with printing speed. Three points in the fracture zone of one tensile test specimen per configuration were measured by DSC analyses. The crystallinity rate varied from 11.4 to 14.3 %, but the standard deviations did not exceed 0.2 %. As observed in Figure 4, a crystallinity rate gradient was revealed between the bottom and the top of the same tensile test specimen. Low standard deviations with three different points in the fracture zone may reveal a constant longitudinal gradient. Models were then significant at a more than 90 % confidence level.

The coefficient of determination for a linear regression with the two mechanical properties as a function of the crystallinity rate was 0.2275 and 0.0269, respectively, indicating a very weak correlation between the mechanical properties and the crystallinity rate. As with the filling rate, the degree of crystallinity was used as a response. The response surfaces have a minimax and then a minimum response from the first to the second design (Figure 9). The same trend was not observed between the mechanical properties and the crystallinity rate as a function of printing speed and extrusion temperature. The printing parameters influence the crystallinity rate, but do not appear to have a major impact on the mechanical properties.

Wittbrodt and Pearce 43 reported that temperature has a relatively linear positive effect on mechanical properties. In fact, this effect appears to be a second-order one. In rheology, an increase in temperature decreases the viscosity and a higher temperature would tend to increase the mechanical properties by developing the adhesion between the filaments 20,24 , due to faster interpenetration kinetics. The contribution of temperature in thermal cycles and thus internal residual stresses 38 and/or thermal degradation due to intra-and inter-molecular exchange, cis-elimination and concerted radical and nonradical reactions above 200 °C44 could have a significant effect and explain the decrease in properties at 215 °C and above in the second Doehlert design. One should note that the addition of clay may contribute to the thermal degradation of the polymer matrix by hydrolysis 45 because of the capacity of this filler to adsorb moisture 40,46 . This phenomenon could be more pronounced for filled materials than unfilled ones if attention is not paid on the drying process before printing. Although temperature has less of an effect than printing speed, the inverse trend of the temperature effect is not explained and could be due to multiple phenomena. Creation of internal residual stresses, ability of materials to be effectively deposited and to weld, modification of the mesostructure, matrix degradation, and crystallisation are all phenomena that can modify the mechanical properties without any being dominant over another. They might express themselves differently according to printing speed and extrusion temperature.

Conclusions and perspectives

Doehlert designs coupled with calorimetry, rheology and tomography were used to better understand the impact of FFF process parameters on tensile properties. The results imply that interpenetration of polymer chains between deposited filaments is not the only major phenomenon that impacts the mechanical properties in terms of thermo-kinetic parameters, especially with composites as interpenetration is reduced by the presence of a filler. Changes in those parameters can both increase and decrease the mechanical properties because of bond formation, abundance of residual stress, crystallisation process, and deposition ability. Because of this complexity of thermo-kinetic interaction, the manufacturing of printed parts should be optimised without dissociation of the temperature and printing speed. The results also showed that rheology under flow modifies the way the filament is deposited, before interpenetration at a zero-shear rate, which in turn impacts the mesostructure. The mechanical properties are then strongly dependent on both linear and non-linear rheology, which implies a certain knowledge of the printed materials behaviour, particularly for composites and the effect on rheology of their filler under flow as viscosity is increased at low shear rate but decreased at high shear rate. Filled materials appeared more suitable for printing at moderate and high printing speed, around 75 mm.s -1 , compared to unfilled materials that are more suitable for low printing speed. Mesostructure optimisation through this DOE might be of interest for architectural control in bone tissue engineering for example.

RSM is recommended for the optimisation of mechanical properties because of the complexity and number of phenomena. The iterative nature of Doehlert designs appeared to be suitable for observing multiple optima for optimisation. This study revealed the efficiency of such method with additive manufacturing, which is recommended for complex system whose the evolution of properties with manufacturing parameters does not correspond to a second-order model like central composite or Box-Behnken designs. Different models are then observed on an entire domain with the use of the Doehlert method, which increases the resolution of the factor effects and broadens the understanding of the results and their physical explanations.

Acknowledgements

The

Figure 1 .

 1 Figure 1. Graphical representations from Makerbot Print Sotfware of a 0° deposition layer (a), a 90° deposition layer (b) and a complete sample (c). Green stripes represent the deposition orientations. Black arrows represent tensile direction.

Figure 2 .

 2 Figure 2. Graphical representation of the two Doehlert designs in coded units with printing speed on the x-axis and extrusion temperature on the y-axis (first design on the left, second on the right). Numbered dots represent the 10 experiments placed according to their respective configuration and the arrow shows the direction of the slip between the first and second Doehlert designs.

Figure 3 .

 3 Figure 3. Response surface for Young's modulus (a, b) and ultimate tensile strength (c, d) as responses for the first (a, c) and second Doehlert theoretical models (b, d) as a function of printing speed and temperature in coded units with isoresponse surface projection.

Figure 4 .

 4 Figure 4. Comparison of DSC thermograms for the two first layers (red) and the two last layers (black) of tensile specimens printed at 215 °C and 75 mm.s -1 (test N° 1).

Figure 5 .

 5 Figure 5. Viscosity of a PLA/layered silicate nanocomposite plotted against shear rate at different temperatures as determined by capillary rheology after Bagley and Rabinovitsch corrections.

Figure 6 .

 6 Figure 6. Images from microCT of printed sample (grip zone) with a cross-sectional view in the z-axis (a, b), x-axis (c, d), and y-axis (e, f) at configuration N° 1 (a, c, e: 75 mm.s -1 at 215 °C) and configuration N° 4 (b, d, f: 15 mm.s -1 at 215 °C). Circles enlighten the observed defects.

Figure 7 .

 7 Figure 7. Images from microCT of fracture zones with a cross-sectional view on the y-axis of tensile specimens printed at 215 °C and (a) 15 mm.s -1 (test N° 4), (b) 45 mm.s -1 (test N° 7), (c) 75 mm.s -1 (test N° 1), and (d) 105 mm.s -1 (test N° 9). Circles enlighten the observed defects.

Figure 8 .

 8 Figure 8. Response surface for (a) the first Doehlert design and (b) the second design with the filling rate as a response as a function of printing speed and the temperature in coded units.

Figure 9 .

 9 Figure 9. Response surface for (a) the first Doehlert design and (b) the second design with the crystallinity rate as a response as a function of printing speed and the temperature in coded units.

Table 1 .

 1 Properties of PLE 005 provided by the supplier.

	Properties	Value
	Specific gravity	1.26
	MFI (190 °C; 2.16 kg)	5-9 g/10 min
	Glass transition temperature	60 °C
	Melting temperature	170 °C-180 °C

The layered silicate was an organically modified montmorillonite (MMT) supplied by Southern Clay Products (Gonzales, TX) under the trade name Cloisite 30B (C30B). The modifier was methyl, tallow, bis-2-hydroxyethyl, and quaternary ammonium with a concentration of 90 mEq per 100 g. The specific gravity of C30B is 1.98, and that of native MMT is 2.86.

Table 3 .

 3 Mean values of Young's modulus and ultimate tensile strength responses as determined by tensile tests for each configuration, with corresponding printing speed and extrusion temperature values, from both Doehlert designs with associated standard deviations.

	Test N°	Printing (mm.s -1 ) speed	Extrusion (°C) temperature	E (MPa)	σM (MPa)
	1	75	215	1397 ± 11	59.4 ± 1.4
	2	60	240	1338 ± 13	56.3 ± 0.6
	3	30	240	1295 ± 41	53.0 ± 2.0
	4	15	215	1345 ± 27	56.7 ± 0.9
	5	30	190	1302 ± 19	56.9 ± 1.7
	6	60	190	1354 ± 19	58.9 ± 1.5
	7	45	215	1249 ± 32	59.0 ± 2.9
	8	90	190	1337 ± 12	58.4 ± 0.4
	9	105	215	1362 ± 18	58.5 ± 0.9
	10	90	240	1357 ± 11	56.3 ± 2.1

Model coefficients after matrix calculations are given in Table

4

. Three-dimensional representations of the response surfaces corresponding to each mathematical model are displayed in Figure

3

. Details of the ANOVA are summarised in the Appendix (Tables

Table 4 .

 4 Coefficients and residue from the theoretical model for both Doehlert designs using Young's modulus and ultimate tensile strength as responses.

			First design		Second design
		Coefficient
			E	σM	E	σM
	Constant	a0	1249.0	59.00	1397.0	59.40
	Printing speed	a1	33.1	1.78	38.0	-0.25
	Temperature	a2	-6.6	-1.88	1.2	-1.36
	Printing speed x temperature	a12	-5.2	0.75	20.8	0.29
	Printing speed x printing speed	a11	122.0	-0.95	-91.5	-0.65
	Temperature x temperature	a22	57.0	-3.32	-36.8	-2.35
	Residue	ε	7.2	0.43	18.5	0.01

Table A1 .

 A1 project PolyFabAdd is co-funded by the European Union. Europe invests in Champagne-Ardenne with the European Regional Development Fund. The authors gratefully acknowledge the Région Champagne-Ardenne (France) [grant number D201501201] and the European Union for their financial support [grant number D201505025]. This project is also included in the MATUR Chair, co-financed by the Région Champagne-Ardenne (France) [grant number D201507798], the European Union [grant number D201507799] and UIMM [Partners convention n°13-2015]. We would also like to thank Mr Philippe Dony and Mr Stéphane Buet for their technical support in this project, and Mrs Annie Leclerc for the proofreading. Analysis of variance of the coefficients from the first Doehlert models using Young's modulus and ultimate tensile strength as responses.

	6. Appendix
	A. ANOVA -First Doehlert design
	Response	Coefficient	Value	Standard deviation	t-value	p-value
		𝑎 0	1249	4.3289	288.526	0.0022
		𝑎 1	33.1333	4.3289	7.6540	0.0827
	𝑬	𝑎 2 𝑎 12	-6.6395 -5.1962	4.3289 4.3289	-1.5338 -1.2003	0.3678 0.4422
		𝑎 11	121.9500	4.3289	28.1712	0.0226
		𝑎 22	57.0167	4.3289	13.1712	0.0482
		𝑎 0	59	0.2907	202.967	0.0031
		𝑎 1	1.7833	0.2907	6.1349	0.1029
	𝝈 𝑴	𝑎 2 𝑎 12	-1.8764 0.7506	0.2907 0.2907	-6.4550 2.5820	0.0978 0.2352
		𝑎 11	-0.9500	0.2907	-3.2681	0.1890
		𝑎 22	-3.3167	0.2907	-11.4097	0.0557

Table A2 .

 A2 Analysis of variance of the theoretical models from the first Doehlert design using Young's modulus and ultimate tensile strength as responses.

		Total	32491.9867	34
		Model	28.6276	5	5.7255	1.9360	0.1188
	𝝈 𝑴	Error	85.7667	29	2.9575
		Total	114.3943	34
		B. ANOVA -Second Doehlert design
	Response	Source of variation	Sum of squares	ddl	Mean square	F-value	p-value
	𝑬	Model Error	13471.5783 19020.4083	5 29	2694.3157 655.8761	4.1080	0.0061

Table B1 .

 B1 Analysis of variance of the coefficients from the second Doehlert models using Young's modulus and ultimate tensile strength as responses.

	Response	Coefficient	Value	Standard deviation	t-value	p-value
		𝑎 0	1397	4.3191	323.446	0.0020
		𝑎 1	38.0000	4.3191	8.7981	0.0720
	𝑬	𝑎 2 𝑎 12	1.1547 20.7846	4.3191 4.3191	0.2673 4.8122	0.8337 0.1304
		𝑎 11	-91.5000	4.3191	-21.1849	0.0300
		𝑎 22	-36.8333	4.3191	-8.5298	0.0743
		𝑎 0	59.4	0.2619	226.772	0.0028
		𝑎 1	-0.2500	0.2619	-0.9544	0.5148
	𝝈 𝑴	𝑎 2 𝑎 12	-1.3568 0.2887	0.2619 0.2619	-5.1798 1.1021	0.1214 0.4691
		𝑎 11	-0.6500	0.2619	-2.4815	0.2439
		𝑎 22	-2.3500	0.2619	-8.9716	0.0707

Table B2 .

 B2 Analysis of variance of the theoretical models from the second Doehlert design using Young's modulus and ultimate tensile strength as responses.

	Response	Source of variation	Sum of squares	ddl	Mean square	F-value	p-value
		Model	10430.5000	5	2086.1000	3.1951	0.0203
	𝑬	Error	18934.5900	29	652.9169
		Total	29365.0900	34
		Model	3.8686	5	1.9737	0.8219	0.5442
	𝝈 𝑴	Error	69.6400	29	2.4014
		Total	79.5086	34


Delineation of cities based on scaling properties of urban patterns: a comparison of three methods

Keywords: urban boundary, urban patterns, scaling laws, fractal, percolation

Identifying urban boundaries involves analysing both the functional and the morphological aspects of urban systems. In this paper, we adopt a purely morphological approach and compare three methods for the morphological delineation of cities. Each method avoids using any predefined quantified threshold (size, distance, built density, etc.) to detect crucial discontinuities in space. The first method identifies Natural Cities by clustering points. The other two are the fractal-based MorphoLim method and the Hierarchical Percolation; both involve transforming the data using a step-by-step dilation process. The three methods are critically compared and illustrated by applications to theoretical urban patterns. We further apply each method to the urban agglomeration of Brussels, the Belgian capital, using different data (building footprints, building centroids and street nodes) and considering two study areas, namely the former province of Brabant and the entire country of Belgium. The results show that it is impossible to draw an unambiguous morphological boundary for an urban agglomeration. Consequently, it is crucial to relate the data used, the size of the study area and the method chosen to the objectives of the delineation.

Introduction

Defining city boundaries is invariably an enthralling challenge for planners and geographers.

Planners need to know where cities start and end in order to govern and manage them efficiently (count populations, levy taxes, organise transport systems, etc.); geographers need to know the city limits in order to map urban growth, characterise land-use changes, compare cities (their areaand population), etc. Sotomayor-Gómez and Samaniego (2020) quote a United Nations World Cities report, which indicates that countries' urban proportions face multiple non-comparable definitions.

As an example, the city of Toronto in Canada may be considered to have between 2.6 and 5.1 million inhabitants, depending on which of the three definitions of city is used to calculate its size.

According to [START_REF] Sotomayor-Gómez | City limits in the age of smartphones and urban scaling[END_REF], such variable figures reflect annual growth rates ranging from 0.9% to 1.8%, which may clearly lead to significant complications for efficient planning efforts. Moreover, in most disciplines, assessing the effect of system boundaries is important to provide certainty and robustness in model findings. This has long been an issue for geographers, from either a theoretical or an empirical perspective as explained in [START_REF] Thomas | City delineation in European applications of LUTI models: review and tests[END_REF].

Even though the "bounding problem" was clearly defined by Haggett (1972) early on, nowadays there is no consensus on the most appropriate way to delineate urban agglomerations either in terms of methods or in terms of criteria or thresholds, and hence results do not converge (see [START_REF] Cottineau | Defining urban clusters to detect agglomeration economies[END_REF][START_REF] Ferreira | Urban settlements delimitation in low-density areas -An application to the municipality of Tomar (Portugal)[END_REF][START_REF] Rozenfeld | The area and population of cities: New insights from a different perspective on cities[END_REF]Thomas et al., 2018, for examples and reviews).

A first way to delineate urban areas is to consider morphological criteria alone based on physical data (buildings, road networks, etc.) (e.g. [START_REF] Medda | Recognition and classification of urban shapes[END_REF][START_REF] Thomas | City delineation in European applications of LUTI models: review and tests[END_REF]. A second way is to consider functional criteria based on variables characterising either places and people living in these places, or exchanges (flows) between places [START_REF] Blondel | Regions and borders of mobile telephony in Belgium and in the Brussels metropolitan zone[END_REF][START_REF] Guérois | Functional Urban Areas (FUA) and European harmonization. A feasibility study from the comparison of two approaches: commuting flows and accessibility isochrones[END_REF][START_REF] Kauffmann | Delineation of City Regions Based on Commuting Interrelations: The Example of Large Cities in Germany[END_REF][START_REF] Sotomayor-Gómez | City limits in the age of smartphones and urban scaling[END_REF]. Here, we have decided to adopt a morphological approach because the criteria used are often easily comparable internationally [START_REF] Fang | Urban agglomeration: an evolving concept of an emerging phenomenon[END_REF]Gisbert et al., 2018;[START_REF] Parr | Spatial Definitions of the City: Four Perspectives[END_REF][START_REF] Weber | An Extended Minimum Spanning Tree method for characterizing local urban patterns[END_REF]. In comparison, functional approaches are less inclusive and capture only part of the urban functioning, depending on the variables considered [START_REF] Thomas | City delineation in European applications of LUTI models: review and tests[END_REF].

Methods used to delineate cities by morphological criteria can be grouped into four categories:

1. Density measures -a density threshold is set, above which a place is considered "urban". This threshold can be applied on discrete maps, where elementary spatial units are administrative units or cells of a grid, or on smoothed maps using a kernel density estimation [START_REF] Borruso | Network Density and the Delimitation of Urban Areas[END_REF]Oliveira et al., 2014;[START_REF] Zhou | Comparative Study of Approaches to Delineating Built-Up Areas Using Road Network Data: Delineating Built-Up Areas Using Road Network Data[END_REF][START_REF] Zhou | Empirical approach to threshold determination for the delineation of built-up areas with road network data[END_REF].

2. Classification of remote sensing images or aerial photographs, so that similar spatial units (pixels) can be grouped together. Usually a criterion of contiguity-generally involving a distance threshold-is added (Abed & Kaysi, 2003;[START_REF] Nicolau | Automated delimitation of urban areas comprising small-sized towns -Comparison of two methodologies applied to mainland Portugal[END_REF][START_REF] Weber | An Extended Minimum Spanning Tree method for characterizing local urban patterns[END_REF].

3. Clustering methods that group elements (usually buildings) according to their proximity or other predefined common characteristics [START_REF] Caruso | Measuring urban forms from inter-building distances: Combining MST graphs with a Local Index of Spatial Association[END_REF]Chaudhry & Mackaness, 2008;[START_REF] Deng | An adaptive spatial clustering algorithm based on delaunay triangulation[END_REF][START_REF] Maisonobe | Delineating urban agglomerations across the world: a dataset for studying the spatial distribution of academic research at city level[END_REF][START_REF] Rozenfeld | Laws of population growth[END_REF][START_REF] Rozenfeld | The area and population of cities: New insights from a different perspective on cities[END_REF][START_REF] Usui | A bottom-up approach for delineating urban areas minimizing the connection cost of built clusters: Comparison with top-down-based densely inhabited districts[END_REF]Wu et al., 2017).

4. Analyses of a statistical distribution across scales -for example, distribution of distances between buildings, street blocks, road intersections or points of interest -and identification of a threshold in this scaling distribution (Arcaute et al., 2016;Jiang, 2013[START_REF] Jiang | Head/tail breaks for visualization of city structure and dynamics[END_REF][START_REF] Jiang | Zipf's law for all the natural cities in the United States: a geospatial perspective[END_REF][START_REF] Jiang | Scaling of geographic space from the perspective of city and field blocks and using volunteered geographic information[END_REF][START_REF] Jiang | The evolution of natural cities from the perspective of location-based social media[END_REF][START_REF] Liu | Delimitating the Natural City with Points of Interests Based on Service Area and Maximum Entropy Method[END_REF][START_REF] Long | Redefining Chinese city system with emerging new data[END_REF][START_REF] Masucci | On the problem of boundaries and scaling for urban street networks[END_REF][START_REF] Tannier | A Fractal Approach to Identifying Urban Boundaries[END_REF][START_REF] Tannier | Defining and characterizing urban boundaries: A fractal analysis of theoretical cities and Belgian cities[END_REF].

The present paper is anchored in the fourth category of methods that avoids using any predefined quantified threshold (size, distance, built density, etc.) to detect crucial discontinuities in space.

Thresholds are identified a posteriori; they are specific to each area under study and cities characterised by similar global densities may exhibit different thresholds. Avoiding ex ante thresholds is especially valuable, as it makes it possible to compare cities of varying sizes and shapes (Le [START_REF] Gléau | Towns of Europe: To Each Country Its Definition[END_REF]. The first method under consideration derives Natural Cities by clustering points (street nodes or building centroids). It combines several methodological processes proposed in a series of five papers [START_REF] Jiang | Zipf's law for all the natural cities in the United States: a geospatial perspective[END_REF], Jiang, 2013, 2015, 2019;[START_REF] Jiang | The evolution of natural cities from the perspective of location-based social media[END_REF]. Two thresholds are successively identified: first, the mean distance between the points of interest, which enables us to create clusters of points; second, the mean area of the clusters of points calculated via the recursive head/tail breaks method applied to the rank-size distribution of the clusters. The Natural Cities identified on this basis are the clusters whose area is larger than the nested means. The other two methods under consideration are the Hierarchical Percolation (Arcaute et al., 2016;[START_REF] Masucci | On the problem of boundaries and scaling for urban street networks[END_REF] and fractal-based MorphoLim methods [START_REF] Tannier | A Fractal Approach to Identifying Urban Boundaries[END_REF][START_REF] Tannier | Defining and characterizing urban boundaries: A fractal analysis of theoretical cities and Belgian cities[END_REF]. Both involve transforming the data (buildings mapped in 2D or road intersections) using a step-by-step dilation. The resulting dilation curves represent the size of the dilation buffer on the x-axis and the change in the number of clusters (in the case of MorphoLim) or the number of road intersections in the biggest cluster (in the case of Hierarchical Percolation) on the y-axis. A crucial threshold is then identified on the dilation curve. This threshold corresponds to a break point in the fractal organisation of the built pattern (MorphoLim) or to the condensation threshold of a logistic function (Hierarchical Percolation).

Few papers broach the comparison of the morphological limits of cities resulting from different methods [START_REF] Borruso | Network Density and the Delimitation of Urban Areas[END_REF][START_REF] Cetinkaya | Proximity-based grouping of buildings in urban blocks: a comparison of four algorithms[END_REF][START_REF] Zhou | Comparative Study of Approaches to Delineating Built-Up Areas Using Road Network Data: Delineating Built-Up Areas Using Road Network Data[END_REF] and none tackle the comparison of the three methods cited above. This is the research objective of this paper. It is well-known that case studies in geography are highly dependent on the size, shape and position of the spatial units used for the analysis (MAUP: Modifiable Areal Unit Problem) [START_REF] Briant | Dots to boxes: do the size and shape of spatial units jeopardize economic geography estimations[END_REF][START_REF] Openshaw | The Modifiable Areal Unit Problem. Concepts and Techniques in Modern Geography[END_REF], as well as the spatial extent of the study area [START_REF] Goodchild | Laws in Geography[END_REF]. Therefore, we seek to determine the impact of these effects on the delineation of morphological agglomerations for each of the three methods.

Data used are of two types: the footprint of buildings (2D) or their centroid (1D), and the road intersections (1D). Such data are easy to collect, freely available in many countries and often easy to handle. This paper is structured as follows: a critical analysis of each method based on the analysis of theoretical urban patterns (Section 2), the application of the three methods to delineate the morphological agglomeration of Brussels (Belgium) (Section 3), and a conclusion discussing the impossibility of defining a single unambiguous morphological limit of a city (Section 4). The theoretical and empirical urban patterns analysed in sections 2 and 3 are freely available in the form of shapefiles (see section entitled "Data and codes availability statement" at the end of this paper).

Description and comparison of three delineation methods based on the analysis of theoretical urban patterns

Natural Cities

The data used to identify Natural Cities are street nodes including intersections and end points [START_REF] Jiang | Zipf's law for all the natural cities in the United States: a geospatial perspective[END_REF]) (Figure 1a), although other types of points could be used, e.g. social media users' locations [START_REF] Jiang | The evolution of natural cities from the perspective of location-based social media[END_REF]. The methodological process adopted to delineate Natural Cities begins with the creation of a Triangular Irregular Network (TIN) as proposed by [START_REF] Jiang | The evolution of natural cities from the perspective of location-based social media[END_REF] (Figure 1b). On this basis, we plot the rank-size distribution of the edges of the TIN according to their length and we verify that this distribution is heavily right-skewed, with a minority of large values in the head and a majority of small values in the tail. If the condition is confirmed, we apply the head-tail division rule set out by Jiang & Liu (2011): "given a variable X, if its values x follow a heavy tailed distribution, then the mean (m) of the values divides all the values into two parts: a high percentage in the tail, and a low percentage in the head" (Figure 1c). Thus, edges shorter than the mean length are selected (Figure 1d). We then transform the sets of contiguous edges into polygons, i.e. morphological clusters. We plot the rank-size distribution of the clusters according to their area. On the condition that this distribution is heavily right-skewed, we apply the head/tail breaks method proposed by Jiang (2013[START_REF] Jiang | Head/tail breaks for visualization of city structure and dynamics[END_REF][START_REF] Jiang | A recursive definition of goodness of space for bridging the concepts of space and place for sustainability[END_REF]. The values of the distribution are separated into two groups: those above the mean and those below the mean; the partitioning for the values above the mean then continues iteratively until the head part of the distribution is no longer heavy-tailed, i.e. more than 40% of the clusters have an area higher than the mean [START_REF] Jiang | A recursive definition of goodness of space for bridging the concepts of space and place for sustainability[END_REF]. By the end of this process, the largest clusters across scales (i.e. the Natural Cities) are identified. In order to select the shortest edges of the TIN, we apply the non-recursive head-tail division rule and not the recursive head/tail breaks method because the rank-size distribution of the edges of the TIN is short-tailed and not long-tailed for the two types of data used and the two study areas considered in Section 3 below. This distinction between short-tailed and long-tailed distributions is pointed out in [START_REF] Jiang | A recursive definition of goodness of space for bridging the concepts of space and place for sustainability[END_REF]: "if the notion of far more smalls than larges recurs on at least twice, the data are long-tailed (...); otherwise, they are short-tailed." In other words, if a distribution is no longer heavy-tailed after the first head/tail break, this means that it is short-tailed. Thus, we use the mean to select the shortest edges of the TIN. For right-skewed data distributions, the median is substantially smaller than the mean. Incidentally, using the median instead of the mean would lead to the creation of smaller Natural Cities. In Figure 2, it can be seen that the delineation of the Natural Cities is more relevant with the median than with the mean. Nevertheless, this observation should be confirmed in the case of real-world urban patterns (see Section 3.5).

The delineation of Natural Cities also crucially depends on the geographical extent of the study area. Expanding the study area may lead to the introduction of supplementary "rural" street nodes being more distant from each other than "urban" street nodes: this mechanically increases the mean length of the edges of the TIN (Figure 3b). Conversely, the introduction of new "urban" street nodes shortens the mean length of the edges of the TIN (Figure 3c).

Finally, concerning the selection of the largest clusters across scales (i.e. the Natural Cities), it is worth noticing that Wang et al. (2015) and [START_REF] Liu | Delimitating the Natural City with Points of Interests Based on Service Area and Maximum Entropy Method[END_REF] propose to use the maximal entropy instead of the head/tail breaks to identify the optimal size of clusters. This alternative is not considered in this paper but could be worth testing further.

MorphoLim

Data used to delineate morphological agglomerations with MorphoLim (https://sourcesup.renater.fr/ www/morpholim/) are building footprints mapped in 2D. A step-by-step dilation is applied to each building on the map; the number of built clusters is counted after each dilation step and the results are presented on a log-log plot, where the x-axis represents the width of the dilation buffer and the y-axis the corresponding number of built clusters (Figure 4a). If a pattern is fractal, then the number of clusters is related to their size by a power-law function corresponding to a linear relationship on a log-log plot and the slope of the curve is the fractal dimension [START_REF] Tannier | A Fractal Approach to Identifying Urban Boundaries[END_REF]. If the pattern is not purely fractal or even multifractal, the dilation curve is not linear and the point at which the dilation curve deviates most from a straight line is a crucial distance threshold. On the curvature function of the dilation curve, which measures how far the curve deviates from a straight line at each point [START_REF] Lowe | Organization of smooth image curves at multiple scales[END_REF], it is possible to identify the point characterised by the maximum curvature (Figure 4b). The maximum curvature of the dilation curve reveals a major spatial discontinuity across scales. The corresponding distance threshold separates two morphological spatial subsets that are distinct in fractal terms: below that threshold, built elements are organised according to the same spatial logic and belong to the same morphological subset. Mapping the urban boundaries then consists in applying to the building map a buffer with a diameter equal to the distance threshold. On this map, the largest built cluster(s) is (are) identified by visual analysis of the rank-size distribution of all built clusters. The largest built cluster(s) is (are) considered to form the urban morphological agglomeration [START_REF] Tannier | Defining and characterizing urban boundaries: A fractal analysis of theoretical cities and Belgian cities[END_REF]. Yet the separation of two morphological spatial subsets that are distinct in fractal terms, which is indicated by the maximum curvature of the dilation curve, does not always correspond to the separation between the city and its periphery, as Figure 5 shows. Figures 5a and5b represent two theoretical built patterns with a densely built-up urban core surrounded by a less dense suburban area that is itself surrounded by an even less dense rural area. The two Figures are really similar except that the urban core in Figure 5a is denser than that in Figure 5b. As a consequence, the maximum curvature in Figure 5a indicates the boundary of the urban core, while in Figure 5b it indicates the boundary of the suburban area. In Figure 5a, the built shape of the suburban area is actually more similar to the built shape of the rural area than to the shape of the urban core. In other words, the strongest morphological difference concerns the urban core with respect to its suburban and rural periphery. Thus, in this case, the delineation of the city excludes the suburban area. In contrast, in Figure 5b, the built shape of the suburban area is closer to that of the urban core.

Consequently, the delineation of the city includes the suburban area.

The delineation of a morphological agglomeration using MorphoLim can depend on the geographical extent of the study area. Figure 6 illustrates this phenomenon. If we look at Figures 6a and6b, the delineation of the morphological agglomeration does not change with the expansion of the study area: the morphological agglomeration comprises the dense urban core and the suburban area surrounding it. This is because the expansion of the study area in Figure 6b with respect to Figure 6a consists in the addition of supplementary buildings in the rural area. Since the spatial distribution of those supplementary buildings conforms with the spatial distribution of buildings located in the rural area of Figure 6a, the morphological difference between the morphological agglomeration and its "rural" periphery is the same in both cases and their delineation is identical.

In contrast, the expanded study area in Figure 6c, compared with Figure 6b, contains a second dense urban core. Thus, the strongest morphological difference is now between, on the one hand, the dense urban cores and, on the other hand, their suburban and rural peripheries. As a result, the suburban area, which surrounds the urban core located on the left of Figure 6c, is no longer included in the morphological agglomeration.

Hierarchical Percolation

Hierarchical Percolation applies to street intersections (Arcaute et al., 2016;[START_REF] Masucci | On the problem of boundaries and scaling for urban street networks[END_REF] (Figure 7). A street network is defined as a planar graph where the street intersections are the vertices and the street segments are the edges. Thus, intersection points are connected only if there is a street connecting them. Two street intersections belong to the same morphological cluster if the distance between them is below a given distance threshold. Progressively increasing the distance threshold via a dilation process enlarges the size of the clusters, until eventually a giant component appears, spanning the entire street network. The size of the largest cluster N max (τ), measured in terms of the number of street intersections, increases with respect to the dilation distance (τ) according to a logistic growth function, where C is the carrying capacity, r the growth rate and τ 0 the inflection point (equation 1). The R code written to calculate the carrying capacity is provided (see section entitled "Data and codes availability statement").

N max ( τ )= C 1+ e -r (τ-τ 0 ) (1)
The city condensation threshold is defined as the threshold where the maximum cluster size N max (τ) intersects the carrying capacity of the logistic function. The largest cluster identified at this condensation threshold forms the city.

Comparison of the methods

All three methods involve the use of fully disaggregated data and follow a bottom-up approach.

Thus, they avoid all biases associated with the MAUP. Nevertheless, their results are influenced by the existence of natural barriers, especially wide rivers or forests, or major transport routes. In the case of Natural Cities, the distances between street nodes separated by any such barrier can be above the mean (Figure 9a). In the case of MorphoLim, a barrier can separate a morphological agglomeration into two clusters. In Figure 9b, the two clusters have almost the same size. Yet if one of these clusters were smaller, it would not appear at the top of the rank-size distribution of built clusters and so would not be identified as a part of the morphological agglomeration. In the case of Hierarchical Percolation, the barriers are reflected as jumps in the logistic function (Figure 9c). To overcome this limitation, [START_REF] Masucci | On the problem of boundaries and scaling for urban street networks[END_REF] have artificially added street intersections in the middle of the bridges on the River Thames in London (UK). Interestingly, city delineation obtained using the three methods can take into account ribbons of buildings or street nodes linking two cities (Figure 10). Finally, each method involves some arbitrariness. As for Natural Cities, the head/tail breaks are defined here as more than 40% of the clusters with an area greater than the mean, but another percentage (20% or 30%) could be applied. As for Hierarchical Percolation, the arbitrariness refers to the visual inspection of the curve representing the growth of the size of the largest cluster in the course of dilation to identify the condensation thresholds. As for MorphoLim, the arbitrariness refers to the choice of the appropriate polynomial to fit the empirical dilation curve before calculating the curvature function.

Real case study: delineation of the morphological agglomeration of

Brussels (capital city of Belgium)

The three methods under consideration have been applied to delineate the morphological agglomeration of the Brussels metropolitan area. While Brussels is undoubtedly the name of a city, its spatial delineation is not clear. From an administrative and political standpoint, Brussels is referred to as Brussels Capital Region (BCR) (Figure 11a). Yet the city extends over a much larger area across the borders of BCR into the Flemish and Walloon regions. The boundaries of this large metropolitan area are not officially or scientifically fixed: indeed, there is no governance at the scale of the metropolitan area and scientists do not concur about its limits (see, for example, [START_REF] Dujardin | Quelles frontières pour Bruxelles ? Une mise à jour[END_REF][START_REF] Boussauw | Relationship between Spatial Proximity and Travel-to-Work Distance: The Effect of the Compact City[END_REF][START_REF] Thomas | Revisiting the extension of the Brussels urban agglomeration: new methods, new data … new results[END_REF][START_REF] Vandermotten | Adéquations et inadéquations du découpage territorial aux contextes politiques et économiques : le cas de la Belgique[END_REF]Van Hecke et al., 2007). Hence the delineation of the morphological agglomeration of Brussels raises not only interesting political issues but also the question of statistical homogeneity throughout the three administrative regions concerned.

Data and study areas

The Hierarchical Percolation and MorphoLim methods necessitate adopting a large study area comprising at least one urban agglomeration and its hinterland (i.e. suburban or rural areas under the influence of the urban core). The Natural Cities method works better with big data, hence a large study area, as the statistical distribution of the data is more likely to be heavy-tailed [START_REF] Jiang | Head/tail breaks for visualization of city structure and dynamics[END_REF]. Additionally, we have shown that all three methods are sensitive to the geographical extent of the study area. For those reasons, we delineate the morphological agglomeration of Brussels within two different study areas: (1) the former province of Brabant [START_REF] Adam | Cartographies des champs d'interaction dans et autour de Bruxelles : navettes, déménagements et appels téléphoniques[END_REF] and (2) the entire country (whole of Belgium). The former province of Brabant (3357 km²) was created in 1815 and is today divided into three entities: BCR (1.2 million inhabitants), the province of Flemish Brabant in the north (1.1 million inhabitants) and the province of Walloon Brabant in the south (0.4 million inhabitants). It includes about ten regional cities (Figure 11b). The entire country of Belgium (30689 km²) has 11.46 million inhabitants (2019).

Data about buildings come from the Land Registry Administration of Belgium (© 2009

Administration Générale de la Documentation Patrimoniale). In this registry, all buildings with a footprint of more than 4 m 2 are mapped in 2D (polygons), regardless of their function (e.g. residential, commercial or industrial). For our study, to avoid noise due to small buildings such as garages or garden sheds, all entities under 12 m², which represent less than 5% of all buildings in the whole of Belgium, were removed from the database. Road network data come from the Geofabrik platform of OpenStreetMap (http://download.geofabrik.de accessed 08/21/2018). OSM data have been chosen because of their large availability across the world. All cycle paths and car parks were removed from the database and all streets with the same name were merged to avoid their over-representation. Within cities, it can happen that street segments are duplicated two or three times in OSM data. Yet a visual inspection of the data has shown that this is rarely the case in Belgium. Street segments have been transformed into street nodes according to the tutorial by [START_REF] Ren | Tutorials for Extracting Street Nodes from OSM Data and Creating Natural Cities from Street Nodes Using ArcGIS Model, A Workshop on Big Data Analytics: Topological and scaling perspective for better understanding and (re-)making sustainable cities[END_REF]. All nodes of degree two were also removed.

Morphological agglomerations of Brussels obtained using the Natural

Cities method

The Natural Cities method is successively applied to street nodes, which include street junctions and end points, and to building centroids (see Appendix 1). The basic condition for applying the Natural Cities method is that the rank-size distribution of the edges of the TIN according to their length is heavily right skewed. This condition is fulfilled both for the province of Brabant and for Belgium.

Deriving Natural Cities from street nodes

When considering both the province of Brabant and the whole of Belgium, a single Natural City forms the morphological agglomeration of Brussels (Figure 12). This Natural City is centred on Brussels Capital Region. Its size is close to that of BCR but its shape differs. The Natural City is slightly larger for the whole of Belgium than for the province of Brabant (148 km² versus 111 km²;

see Table 1).

Deriving Natural Cities from building centroids

When considering the province of Brabant, four small Natural Cities make up the morphological agglomeration of Brussels. Three are located within Brussels Capital Region and the fourth covers the regional cities of Braine-l'Alleud and Waterloo (see Figure 11b and Figure 13a). When considering the whole of Belgium, only one Natural City located within Brussels Capital Region makes up the morphological agglomeration of Brussels (Figure 13b).

Morphological agglomerations of Brussels obtained using MorphoLim

Data used here are buildings mapped in 2D. The analysis was first performed on the province of Brabant. Figure 14a shows that the dilation curve is quite close to a straight line. The maximum curvature value is low (-0.087). The spatial distribution of buildings in the cluster at the top of the rank-size distribution (Figure 14b) therefore differs little from the spatial distribution of buildings in the rest of the province of Brabant. This largest cluster covers a major part of the study area (Figure 15a) and incorporates all regional cities in Figure 11b apart from Tienen and Nivelles.

Indeed, Nivelles is quite confined within itself and Tienen is connected by a few built ribbons to some nearby villages but not to the other regional cities. The results are markedly different when the study area is the whole of Belgium: the dilation curve exhibits a stronger main curvature (equal to -0.11) (Figure 14c); the cluster at the top of the ranksize distribution is not so large (Figure 14d); the morphological agglomeration spreads towards the south (Walloon Brabant) and towards the north-west (Flemish Brabant) beyond the border of the province where it stops at the edge of the city of Ghent (Figure 15b).

Morphological agglomerations of Brussels obtained using Hierarchical

Percolation

Hierarchical Percolation applies to road intersections (without the street end points unlike the Natural Cities method). Adaptation of the method enables us to apply it to building centroids too [START_REF] Behnisch | Settlement percolation: A study of building connectivity and poles of inaccessibility[END_REF].

Hierarchical Percolation applied to road intersections

In the case of the province of Brabant, the first condensation threshold appears at 140 m. The morphological agglomeration of Brussels extends to the close periphery of Brussels Capital Region and a little to the south to Braine-l'Alleud (Figure 16a). With the second threshold (170 m), the agglomeration extends beyond the administrative limits of BCR in all directions and encompasses Braine-l'Alleud in the south. The results obtained for Belgium are very similar (Figure 16b): a first percolation threshold appears at 140 m and the morphological agglomeration of Brussels extends to the close periphery of BCR and to Braine-l'Alleud.

Turning now to the largest condensation thresholds, the morphological agglomeration of Brussels is no longer identifiable. In the case of the province of Brabant (threshold: 220 m), the resulting morphological agglomeration covers almost the entire study area. In the case of Belgium (threshold: 440 m), the resulting morphological agglomeration covers almost the entire country. x-axis: width of the dilation buffer; y-axis: number of road intersections in the largest cluster; τ: condensation threshold.

Hierarchical Percolation applied to building centroids

Similarly to the above-mentioned results, several percolation thresholds emerge at different dilation steps. For the province of Brabant, three thresholds are identified (Figure 17a). If we consider the first one (30 m), the agglomeration of Brussels is almost entirely limited to BCR and only extends beyond it to the south where it encompasses Braine-l'Alleud. Considering the second threshold (90 m), the agglomeration of Brussels is much larger and covers almost all of the western half of the study area. With the third threshold (100 m), the two largest cities in the study area (Brussels and Leuven) merge and the morphological agglomeration covers all regional cities apart from Tienen and Nivelles.

With the extension of the study area to the whole of Belgium, the first percolation threshold appears at 60 m (Figure 17b). At this stage, Brussels covers the western half of the study area and encompasses all regional cities apart from Leuven, Tienen and Nivelles. With the second threshold (70 m), the agglomeration of Brussels encompasses the city of Ghent but still not Leuven. Finally, with the next threshold (160 m), all of Flanders is attached to the agglomeration of Brussels because of the existence of built ribbons that link the Flemish cities.

Discussion

By using three methods for delineating the morphological agglomeration of Brussels, we obtain a series of morphological limits that differ markedly. With Hierarchical Percolation, MorphoLim and Natural Cities applied to street nodes, the morphological agglomeration is made up of a single morphological cluster. In contrast, the application of the Natural Cities method to building centroids results in the identification of several small Natural Cities. In general, all Natural Cities are quite small. It would therefore be irrelevant to use the median instead of the mean, as considered in Section 2.1.

With MorphoLim, the morphological agglomeration of Brussels is smaller if the study area is Belgium (Figure 15b) than if it is the province of Brabant (Figure 15a); see Table 1. This is because a crucial change in distances separating buildings occurs at a lower distance threshold when considering all of Belgium. The opposite is found for the Hierarchical Percolation method applied to building centroids when the first percolation threshold is considered (see Figures 17a andb, and Table 1). Yet in the case of street intersections, the percolation thresholds and the area of the morphological agglomerations are the same whichever study area is under consideration (province of Brabant or whole of Belgium) when the first percolation threshold is considered (see Figures 16a andb, and Table 1). For Natural Cities, the effect of the expansion of the study area also depends on the data used. With the street intersections, we observe an increase in the mean length of the edges of the TIN and a decrease in the minimal area of the clusters. With the building centroids, we again observe an increase in the mean length of the edges of the TIN, but this time an increase in the minimal area of the clusters (see Figures 12 and13, and Table 1).

Among all delineations of the agglomeration of Brussels, some are undoubtedly irrelevant. First, Natural Cities derived from building centroids are much too small (Figure 13). They delineate large neighbourhoods within BCR as well as the morphological agglomeration of two regional cities (Braine-l'Alleud and Waterloo), but not one inclusive morphological agglomeration. Second, the boundaries of the morphological agglomeration are sometimes artificially truncated by the border of the province of Brabant when this study area is chosen. This is the case with the delineation obtained using MorphoLim (Figure 15a) as well as the application of the largest percolation thresholds of Hierarchical Percolation (Figure 16a -threshold 220 m, Figure 17a -threshold 100 m). In those cases, the province of Brabant is too small with respect to the chosen method or parameters and a larger study area (e.g. the whole of Belgium) would be required. Third, when the study area is the whole of Belgium, the largest percolation thresholds do not enable us to delineate the morphological agglomeration of Brussels because the morphological cluster being identified covers either all of Flanders (for building centroids) or almost all of Belgium (for street intersections).

Finally, nine delineations of the morphological agglomeration of Brussels appear to be geographically relevant: the Natural Cities identified on the basis of street nodes for the two study areas (province of Brabant and Belgium), the delineation obtained using MorphoLim applied to the whole of Belgium, the delineations obtained using the Hierarchical Percolation method applied to road intersections (province of Brabant -thresholds 140 m and 170 m, and Belgium -threshold Finally, two opposing views of the settlement pattern of Belgium are supported by the results obtained. On the one hand, the limit of the morphological agglomeration of Brussels remains in the close periphery of BCR or spreads a little into the Walloon region to Braine-l'Alleud (six relevant delineations of Brussels out of nine). This confirms the idea that Belgian cities are still separate entities [START_REF] Van Criekingen | Brussels: polycentricity as 'images on the map[END_REF]. In contrast, in three delineations out of nine, the morphological agglomeration of Brussels extends clearly over the borders of BCR into the Flemish and Walloon regions, or spills over the border of the province of Brabant and reaches the city of Ghent in the north-west. Moreover, the morphological agglomeration even covers either all of Flanders or almost all of Belgium, with the largest percolation thresholds. This supports the idea that Belgium is one large conurbation [START_REF] Vandermotten | European Planning and the Polycentric Consensus: Wishful Thinking?[END_REF] in opposition to [START_REF] Van Criekingen | Brussels: polycentricity as 'images on the map[END_REF]. In actual fact, Belgian cities are separate entities but they have strong morphological ties that take the form of built ribbons, especially in Flanders.

Conclusion

Conventional methods that integrate a fixed threshold set ex ante have the advantage of being easy to apply and their results are easy to handle. Yet such methods overlook the specific features of the patterns under consideration that appear across scales and that are revealed by the three methods compared in this paper. Each of these methods has a specific way of detecting the urban boundary -the nested mean value in a head/tail breaks distribution for Natural Cities, a maximum deviation from a hyperbolic distribution for MorphoLim and a condensation threshold in a logistic growth function for Hierarchical Percolation. Cities with the same built density or the same density of street nodes may display different dilation thresholds (MorphoLim) or percolation thresholds (Hierarchical Percolation), or even mean edge lengths (Natural Cities). All three methods involve the use of fully disaggregated data and follow a bottom-up approach. Consequently they avoid all biases of the MAUP. Nevertheless, they are all sensitive to the geographical extent of the study area.

The three methods were applied to two study areas (the province of Brabant and the whole of Belgium) with two types of data (buildings and street intersections). The resulting delineations of the morphological agglomeration of Brussels differ markedly. Accordingly, it is impossible to draw a single unambiguous morphological boundary for an urban agglomeration in general, and for Brussels in particular. This is consistent with the observations by Chen (2016): "Urban form has no characteristic scale, and thus an urban boundary cannot be identified exactly" and by Jiang (2018): "(...) things are not measurable, or measurement depends on the measuring scales because of the fractal nature of geographic phenomena (e.g. [START_REF] Goodchild | The Fractal Nature of Geographic Phenomena[END_REF][START_REF] Batty | Fractal Cities: A geometry of form and function[END_REF][START_REF] Chen | Derivation of the functional relations between fractal dimension of and shape indices of urban form[END_REF])". Indeed, if a city is a boundless fractal, none of the three methods can enable its boundary to be identified. As suggested by [START_REF] Chen | Defining urban boundaries by characteristic scales[END_REF], "(...) geographical phenomena fall into two groups: one is with characteristic scale (scaleful group), and the other, without characteristic scale (scale-free group). The former can be described with characteristic length such as average value, standard deviation, and eigenvalue, while the latter should be described with scaling exponent such as fractal dimension." Nevertheless, the three methods make it possible to delineate a morphological agglomeration of Brussels because local deviations from scale invariance actually exist [START_REF] Chen | Defining urban boundaries by characteristic scales[END_REF][START_REF] Sémécurbe | Applying two fractal methods to characterise the local and global deviations from scale-invariance of built patterns throughout mainland France[END_REF]. As stated by [START_REF] Jiang | Geographic space as a living structure for predicting human activities using big data[END_REF], geography has three fundamental issues regarding geographical space: how it looks, how it works, and what it ought to be. In the present paper, we compare three methods with a view to seeing how they describe the built-up realities on the basis of the example of the morphological agglomeration of Brussels. However, we are not comparing these methods with a view to finding out how a morphological agglomeration should be, nor are we trying to determine whether or not results comply with norms or ideas of goodness or optimality, as done for instance in Chen (2016) and [START_REF] Jiang | A recursive definition of goodness of space for bridging the concepts of space and place for sustainability[END_REF]. With the example of Belgium, we here compare methodologically and empirically three methods for delineating morphological agglomerations; Belgium is sufficiently diverse in terms of urban textures for trusting the present results. Comparing with case studies in other countries/continents opens a new avenue for further research in human and urban geography.

Our study has shown that street nodes and building centroids do not convey the same information about the spatial organization of built patterns. This is very clear in the case of the Natural Cities method, as the Natural Cities delineated on street nodes are much larger than the Natural Cities delineated on building centroids. Distances that separate street intersections are actually above the mean in between individual cities, whereas distances that separate building centroids are above the mean everywhere except in a few very dense intra-urban built clusters. Moreover, the nature of the data used exerts a joint effect with the size of the study area on city delineations. This joint effect varies with the method under consideration. With this in mind, the crucial interest of the Natural Cities method applied to street nodes is that it enables clear identification of the morphological limits of individual cities, in this case Brussels but it could be other cities as well. Considering the Belgian case, changing the study area has a minor impact on the resulting delineations. Indeed, expanding the study area expands Brussels' Natural City just a little bit. In contrast, the MorphoLim and Hierarchical Percolation methods applied to the building data for the whole of Belgium enable the identification of an extensive morphological agglomeration of Brussels that includes many regional cities, among which Brussels, linked by built ribbons. However, the application of those two methods to building data considering only the province of Brabant is not relevant because the boundaries of the morphological agglomeration are truncated by the border of the study area.

Finally, the Hierarchical Percolation method applied to street intersections is interesting, as it enables the identification of three percolation thresholds, and thus three nested morphological agglomerations of growing size that show how the peripheries of Brussels are more or less closely linked to the core agglomeration from a morphological point of view. Conversely, this method does not provide a single clear limit of the morphological agglomeration of Brussels.

To conclude, geographers and planners should be more than ever conscious of the importance of their choices about the data used, the size of the study area and the delineation method. They all have to be carefully selected with respect to both the objectives of each delineation and the specific features of the area under study.

Figure 1 .

 1 Figure 1. Methodological process adopted to delineate Natural Cities: (a) identification of the street nodes; (b) creation of a Triangular Irregular Network (TIN); (c) rank-size distribution of the edges of the TIN according to their length; (d) selection of the edges shorter than the mean; (e) transformation of the selected edges into polygons; (f) selection of the largest clusters (i.e. the Natural Cities).

Figure 2 .

 2 Figure 2. A road network (a) and two variants of the Natural Cities method: delineation based on the mean length of the edges of the TIN (b); delineation based on the median length of the edges of the TIN (c).

Figure 3 .

 3 Figure 3. Effect of the geographical extent of the study area on the delineation of Natural Cities.

Figure 4 .

 4 Figure 4. MorphoLim: identification of the dilation threshold marking the limit of the morphological agglomeration (blue disks).

Figure 5 .

 5 Figure 5. Analysis of two theoretical built patterns using MorphoLim. The urban core of pattern (a) is denser than the core of pattern (b).

Figure 6 .

 6 Figure 6. Effect of the geographical extent of the study area on the identification of city boundaries using MorphoLim.

Figure 7 .

 7 Figure 7. Hierarchical Percolation method applied to a theoretical street network. (a) identification of street intersections; (b) fitted logistic function (x-axis: size of the dilation buffer; y-axis: number of street intersections in the biggest cluster); (c) city delineation.

Figure 8 .

 8 Figure 8. A polycentric street pattern (two centres) (a), and a double condensation threshold in the logistic growth function of the size of the largest cluster (b).

  All three methods are also sensitive to the geographical extent of the study area. Yet the underlying mechanism is different. For MorphoLim and Hierarchical Percolation, the extent of the study area does not change the dilation threshold or the percolation threshold unless a new built structure appears with the expansion of the study area. In contrast, city delineations obtained using the Natural Cities method vary with respect to the expansion of the study area even if no new built structure is included.In the case of Natural Cities, the distances separating street nodes are measured on a Triangular Irregular Network (TIN). In the case of Hierarchical Percolation and MorphoLim, distances separating buildings are Euclidean distances. City delineation with Natural Cities results from direct analysis of the distribution of distances that separate street nodes. With MorphoLim and Hierarchical Percolation, the distribution of distances is not directly analysed: city delineation appears with a jump in the size of the largest cluster of street intersections in the course of a dilation process (Hierarchical Percolation) and a maximum deviation from scale invariance of the number of built clusters in the course of dilation (MorphoLim). Although the criteria (measurement and analysis of distances) are different between the three methods, they all look for heterogeneity, i.e. spatial differentiation. Consequently, if buildings or street intersections are uniformly distributed in space, i.e. if distances separating buildings do not vary, no morphological agglomeration can be identified, whatever the method chosen. Besides, with Hierarchical Percolation and MorphoLim, no morphological agglomeration can be detected if local variations in distances separating the elements (buildings or street intersections) are always the same across scales.

Figure 9 .

 9 Figure 9. Morphological agglomeration within a theoretical urban pattern crossed by a large non built area: (a) Natural Cities; (b) MorphoLim; (c) Hierarchical Percolation.

Figure 10 .

 10 Figure 10. Morphological agglomeration of a theoretical urban pattern where a built ribbon links two urban centres by means of (a) Natural Cities, (b) MorphoLim, and (c) Hierarchical Percolation

Figure 11 .

 11 Figure 11. Belgium: location of the former province of Brabant around the Brussels Capital Region (a). Smoothed densities of building centroids (b) and road intersections (c) in the province of Brabant. Note: the bandwidths correspond to the optimal values proposed by default in the ArcGis software application.

Figure 12 .

 12 Figure 12. Morphological agglomerations of Brussels obtained using Natural Cities computed on street nodes (a) for the province of Brabant and (b) for the whole of Belgium.

Figure 13 .

 13 Figure 13. Morphological agglomeration of Brussels obtained using Natural Cities on the basis of building centroids: (a) for the province of Brabant and (b) for the whole of Belgium.

Figure 14 .

 14 Figure 14. Dilation curve and rank-size distribution of built clusters obtained using MorphoLim for the province of Brabant (a & b) and for the whole of Belgium (c & d).

Figure 15 :

 15 Figure 15: Morphological agglomerations of Brussels obtained using MorphoLim; (a) for the province of Brabant; (b) for the whole of Belgium.

Figure 16 .

 16 Figure 16. Logistic functions and morphological agglomerations of Brussels obtained using Hierarchical Percolation on the basis of road intersections; (a) for the province of Brabant; (b) for the whole of Belgium.

Figure 17 .

 17 Figure 17. Logistic functions and morphological agglomerations of Brussels obtained using Hierarchical Percolation on the basis of building centroids; (a) for the province of Brabant; (b) for Belgium. x-axis: width

  Figure 17b -threshold 70 m).

Table 1 .

 1 Summary of results obtained using the three delineation methods. (*) Mean length of the TIN;( + ) distance at which the dilation curve deviates most from a straight line; ( # ) percolation threshold.

	Method	Data	Study area	Threshold (in m)	Area of Brussels' agglomeration (in km 2 )	Perimeter of Brussels' agglomeration (in km 2 )
	Natural Cities	Road	Province of Brabant	212 (*)	110	169
		intersections	Whole of Belgium	250 (*)	148	204
		Building	Province of Brabant	48 (*)	67	329
		centroids	Whole of Belgium	54 (*)	46	187
	MorphoLim	Building	Province of Brabant	181 ( + )	653	5299
		footprints	Whole of Belgium	121 ( + )	502	5186
	Hierarchical	Road	Province of Brabant	140 ( # )	271	1038
	percolation	intersections		170 ( # )	395	1262
				220 ( # )	1395	4165
			Whole of Belgium	140 ( # )	269	1030
				440 ( # )	21580	20876
		Building	Province of Brabant	30 ( # )	128	1261
		centroids		90 ( # )	710	3053
				100 ( # )	1176	5093
			Whole of Belgium	60 ( # )	395	2073
				70 ( # )	1158	6818
				160 ( # )	8779	25463

Data and codes availability statement

The data and codes that support the findings of this study are available on figshare.com: https://doi.org/10.6084/m9.figshare.11405598

Appendix 1: Head/tail breaks method applied to the rank-size distribution of morphological clusters, i.e. polygons formed by the continuous edges of the TIN whose length is shorter than the mean edge length. Grey line in the Table: end of the head/tail breaks process. After this step, the percentage of clusters in the head of the distribution is equal to 1 or higher than ~40%.

Street nodes of the province of Brabant

Step