Organization Patterns of Complex River Networks in Chile: A Fractal Morphology

Abstract

River networks are spatially complex systems difficult to describe by using simple morphological indices. To this concern, fractal theory arises as an interesting tool for quantifying such complexity. In this case of study, we have estimated for the first time the fractal dimension of Chilean networks distributed across the country, analysed at two different scales. These networks insert into variable environments, not only from a climatic and hydrological point of view, but also from a morphological point of view. We investigate to which extent the fractal dimension is able to describe the apparent disorganized character of landscape, by applying two methods. Striking patterns of organization related to Horton ratios and the fractal dimension are reported and discussed. This last parameter depends on the scale of the network, showing interesting groupings by tectonic and climatological factors. Our results suggest that under restricted conditions, the fractal dimension could help to capture the intricate morphology of Chilean networks and its links with the hydrological, climatic, and tectonic conditions present across the country.

1. Introduction

Fractals are complex objects whose geometrical structure remains invariant regardless of the observation scale [1]. These objects can be found in different areas such as medicine, physics, mathematics, geology, biology, and particularly, in geomorphology. In this last area, many natural systems can be described by using fractal morphometry, for example, the shape of relief, borderlines of lakes, coasts and rivers and particularly, the topology of streams networks [1,2,3,4]. In contrast to some classical metrics introduced for basins [5], the fractal dimension could help to get deep insight about the complex structure of river systems, providing the topological dimension where the structural invariability of the network is observed [1,6]. To determine to which extent such dimension can describe these structures, the collection, and measuring of real river networks remains a fundamental task. In this context, Chilean river networks arise as an interesting case of study of fractal analysis providing a rich scenario to explore the variability of the fractal dimension and its possible links with the different landscapes, climate, and tectonics characteristics of the country well described by Subercaseux [7]. The networks analyzed in this paper were extracted at two different observation scales, emphasizing the possible links with morphometric properties such as the Horton ratios, the drainage density, the geomorphic indices and the geomorphological processes influencing the shape of the landscape (e.g., tectonics, lithology, and climate, among others).

Horton [8] pioneered on this topic introducing a quantitative description of river networks, considering a hierarchical ordering of the streams. The ordering of a river network is a process of determining the level of branching of these systems [9]. Strahler [10] provided precise rules to determine such ordering for a given river system. In this context, Horton defined the following set of ratios representing the streams development [11]:

$$R_{B,\omega }=\frac{N_{\omega -1}}{N_{\omega }}{R}_{L,\omega }=\frac{L_{\omega }}{L_{\omega -1}}{R}_{A,\omega }=\frac{A_{\omega }}{A_{\omega -1}},$$

(1)

where $\omega =2,\cdots ,\Omega$ is the order of the streams. First order rivers ($\omega =1$) are the outermost tributaries of the system, meanwhile higher order segments ($\omega \to \Omega$) correspond to mainstreams of the network [10]. The variables $L_{\omega }$, $N_{\omega }$ and $A_{\omega }$ correspond to the length and number of drains of the sub-catchments of area $A_{\omega }$ and the parameters $R_{A,\omega },R_{B,\omega }$ and $R_{L,\omega }$ denote the area, bifurcation and length ratios of a river system, respectively. In particular, the bifurcation ratio ($R_B$) describes the branching pattern of a drainage network and is defined as ratio between the total numbers of stream segments of order $\omega -1$ to that of the next higher order $\omega$ [12]. An interesting observation from Horton is that for high enough dense networks, $R_{A,\omega },R_{B,\omega }$ and $R_{L,\omega }$ show almost convergent values, a feature of a self-similar (or fractal) behaviour [13,14]. By applying a recursive argument the relationships shown by Equation (1) can also be written as follows:

$$N_{\omega }={\left(R_{B,\omega }\right)}^{\omega -\Omega }{L}_{\omega }=L_1{\left(R_{L,\omega }\right)}^{\omega -1}{A}_{\omega }=A_1{\left(R_{A,\omega }\right)}^{\omega -1},$$

(2)

From these relationships, the average ratios $R_A,R_B,R_L$ can be obtained by computing the mean slope of the curves $log\left(N_{\omega }\right),log\left({\overline{L}}_{\omega }\right)$ and $log\left({\overline{A}}_{\omega }\right)$ versus $\omega$, respectively, where ${\overline{L}}_{\omega }=L_{\omega }/N_{\omega }$ and ${\overline{A}}_{\omega }=A_{\omega }/N_{\omega }$. According to Kim [15] and Beer and Borgas [16], Horton ratios can be seen as a kind of mathematical transformation, such that the structure of channel networks is transformed into deterministic structures where such ratios are valid, limiting their applicability in real networks showing some degree of departure from deterministic structures. Following this idea other studies are concerned about some interpretative restrictions of these ratios [17]. Despite these cautions, these ratios can still be considered as characteristic parameters of a river network.

Inspired on Hack’s law [18], Mandelbrot explored the connection between fractals and the meandering patterns of rivers suggesting the law $\mathit{l}\propto A^{d/2}$, where A is the drainage area, $\mathit{l}$ the length of the mainstream and $d\approx 1.1$ its meandering fractal dimension. Based on this relation and assuming the drainage density constancy across the network, Feder [6] proposed the law $d=2log\left(R_L\right)/log\left(R_B\right)$. Rosso et al. [19] extended this result suggesting the law $d=max\{1,2log\left(R_L\right)/log\left(R_A\right)\}$, in agreement with the topological minimum path dimension deduced by Liu [20].

Fractal Dimension of River Network: A Fractal Tree Approach

One of the problems of previous formulations is the need of connecting them with the fractal dimension of the entire network. La Barbera and Rosso [2] derived a theoretical law for this situation, assuming that Horton parameters holds through the whole network across different scales:

$$D_1=\frac{log\left(R_B\right)}{log\left(R_L\right)}if{R}_B>R_L\&D_1=1if{R}_B<R_L,$$

(3)

According to La Barbera and Rosso [2], $D_1$ leads to values strictly comprised between 1.5 and 2.0 (1.67 in average). The authors argued about the impossibility of reaching values close to 2, car river networks show decreasing drainage densities for increasing contributing areas. Tarboton et al. [21] argues that La Barbera and Rosso [2] assumes that mainstreams are topologically similar to objects of dimension 1. However, many natural streams show visible patterns of deformation (e.g., meanders, avulsions), for which $d\ne 1$. For such condition, Tarboton et al. [21] proposed a modification of the previous law, $D_2=dlog\left(R_B\right)/log\left(R_L\right)$ where $d\approx 1.14$, limiting $D_2$ to values strictly lower than 2 [3,21,22]. This upper bound is coherent with some observations made at larger scales where it is reasonable to assume that streams drain each point of the basin as pointed by [23]. La Barbera and Rosso [24] refuted the conclusion of Tarboton proposing the modification $D_3=\beta log\left(R_B\right)/log\left(R_L\right)$ where $\beta =1/(2-d)$. Liu [20] worked theoretically with infinite dense networks ($\Omega \to \infty$) at both meso and microscale, proposing different dimensions describing the structure of individual streams, one of them related to fractal dimension of the whole network. The next expression was also derived by [21], and writes as follows:

$$D_4=2\frac{log\left(R_B\right)}{log\left(R_A\right)},$$

(4)

Although $D_1,\cdots ,D_4$ are a practical approach to determine the fractal dimension of river networks, they require a huge amount of geographic information. On the other side, there are some limitations of these methods that deserves to be considered. One of them is related with the self-similarity hypothesis behind a drainage geometry. According to Mandelbrot [25], self-similarity is the property shown by some objects formed by portions that can be considered a reduced-scale image of the whole. This property has been objectively refuted by some authors [17,26], leading to other approaches based on the concept of self-affinity. Self-affinity, by contrast, defines an object formed by portions that do not preserve the shape of the whole. Such portions are scaled by different amounts along the relevant directions of the object. Self-affine pattern of river networks was first underlined by Nikora [27,28,29] and it is frequently related with morphological anisotropy of landscapes [23,30]. According to [31], such anisotropy is issue from the combination of different tectonic processes that constrains the diffusion of the streams. Despite these limitations, the estimators $D_1$ to $D_4$ have proven to give further insight about the fractal structure of a river network, compared to the typical morphometric indices used for watersheds [15,32,33,34,35,36], for example, the form factor, the basin circularity or the elongation index defined in Section 2.2.

2. Methodology

2.1. Estimation of the Fractal Dimension

The determination of the fractal dimension in this study was conducted by following two different approaches. The first method is supported by the laws based on Horton metrics (Equations (2) and (3)) and the second approach, is based on the box-counting method applied through the software Fractalyse. This software was developed at THEMA Laboratory and it can be used to obtain an empirical estimation of the fractal patterns of urban and natural networks [37,38]. In our case, this software was fed with Landsat-5 satellite images obtained from NASA’s platform site. Every image was analysed by using image processing tools in Matlab, allowing to extract the planar structure of each network. The key point of the method is supported by the concept of topological recovering of a surface, that is, a given drainage pattern could be covered by a finite number of $N\left(s\right)$ squared-boxes of side s. The box-size s can be reduced step by step by following an arithmetic or a geometric progression, increasing the number of boxes covering the figure. According to Rodriguez-Iturbe and Rinaldo [39], in this method the fractal dimension of the network ($D_F$) can be calculated from the next relationship:

$$D_F=\underset{s\to 0}{lim}\frac{log\left(N\right(s\left)\right)}{log\left(\frac{1}{s}\right)},$$

(5)

This dimension is typically estimated from the slope of the linear part of the curve $log\left(N\right(s\left)\right)$ vs. $log\left(s\right)$. Such linearity is defined by a range of box-sizes, let us say $[s_{min},s_{max}]$. Beauvais and Montgomery [40] warned about the calculation of $D_F$, proposing a method to properly chose this range. The proposed methodology avoids finite size effects [41], leading to fractal dimensions strictly comprised in the range 1$<D_F<$2. In Fractalyse we can easily chose these points. Given the variable map scales of our images, we have defined $s_{min}$ ranging from 29 m to 360 m and $s_{max}$ from 2 km to 37 km. The minimum box-size ($s_{min}$) agrees with the minimum resolution of every image. Under these considerations, the values of $D_F$ were obtained within a 95% confidence band and its associated mean error ${\Delta }_F$.

2.2. Estimation of Shape Indices

In order to compare the morphological characteristics of each network, we have estimated some dimensionless geomorphic parameters used as rough descriptors of a drainage basin [5], that is, the form factor F, the basin circularity C and the basin elongation E, defined as follows [42]:

$$F=\frac{A}{L^2},C=\frac{4\pi A}{P^2},E=\frac{2\sqrt{A}}{\sqrt{\pi }L},$$

(6)

where L is a characteristic length of the watershed (e.g., basin’s runoff distance) and P its perimeter. According to Horton [8], the factor F indicates the flow intensity of a basin. The value $F=\frac{\pi }{4}\approx 0.79$ represents a perfectly circular basin and $F=1$ a squared-shape basin. High F values describe large peak flows of short duration, whereas elongated basins with low values of F show lower peak flows of longer duration.

By contrast, the basin circularity C is mainly controlled by geology and structure, slope, climate and stream properties [43]. Values of C from 0.4 to 0.7 corresponds to elongated and very permeable homogeneous geologic materials [44]. A high value of C is a sign of a landscape in late maturity stage. The basin elongation E reads as the ratio of the diameter of a circle of the same area as the basin to the maximum basin length [43]. Then, the shape of basins can be also sorted with this index, i.e., circular ($0.9\le E\le 1.0$), oval ($0.8\le E\le 0.9$), less elongated ($0.7\le E\le 0.8$), elongated ($0.5\le E\le 0.7$), and more elongated ($E<0.5$). When $E\approx 1$ basins shows less structural effects and when E decreases from 0.8 to 0.6 the basins are dominated by steep gradients and high elevations [10]. According to Bull [45], the elongation index can be also used to describe possible tectonic activity, ranging from high ($E<0.5$) to slight activity ($0.51\le E\le 0.60$).

The drainage density $\rho =Z/A$ of the network is another parameter provided by GRASS-GIS, where Z is the total length of the streams. The drainage density provides a quantitative measure of the average length of streams for the whole basin [8,10]. High values of $\rho$ is a sign of weak or impermeable subsurface material with sparse vegetation and mountainous relief. By contrast, low values of this parameter can be usually associated with regions of high erosion resistance, highly permeable, and dominated by low reliefs leading to coarse drainage textures [10].

2.3. Definition of the Region of Study

In this study, we have analyzed 23 large-basins, located between the latitudes 17${}^{\circ }{30}^{\prime }$ S and ${56}^{\circ }30{}^{\prime }$ S, covering a total area of 363,354 km${}^2$. The area of each network is larger than 6219 km${}^2$ and they were delimited by following the guidelines proposed by Direccion General de Aguas (DGA), a Chilean governmental agency focused on water resources management. Figure 1a shows the location of these units from north to south. This set of networks are representative of the tectonic, geographic and climatic diversity of the country. The geographical analysis of every unit was conducted on the software GRASS-GIS. This software provides the Horton ratios $R_A,R_B,R_L$, the hierarchical order $\Omega$, catchment and sub-catchments drainage areas A, the average slope $i_m$ and the mean elevation $H_m$ of each network. The order $\Omega$ was determined by following the criteria proposed by Strahler [46]. This information is presented in

Table 1. Morphometric parameters of large basins sorted from north to south and by tectonics influence. (NPN: Nazca plate-north, NPFS: Nazca plate-flat slab, NPS: Nazca plate-south, AP: Antarctic plate). Here A is the drainage area, $i_m$ the mean slope, $H_m$ the mean elevation of the unit and $\Omega$ the Strahler’s maximum order of each network (${}^{*}$ Basins chosen for analysis at sub-basin level).

Tectonics	Basin	A (km2)	${\mathit{i}}_{\mathit{m}}$ (%)	${\mathit{H}}_{\mathit{m}}$ (m)	$\mathbf{\Omega }$	F	C	E	$\mathit{\rho }$ (km−1)
NPN	Loa${}^{*}$ (LO)	51,056	13	2401	9	1.88	0.40	0.27	1.86
	Caracoles (QC)	32,537	9	1947	10	0.84	0.43	0.33	2.15
	Salado (SA)	16,826	20	3086	9	0.28	0.33	0.19	1.55
	Average					1.00	0.38	0.26	1.85
NPFS	Copiapo (CO)	18,608	33	2707	8	0.49	0.34	0.24	1.16
	Huasco (HU)	9759	43	2738	7	0.40	0.31	0.25	0.98
	Elqui${}^{*}$ (EQ)	9484	46	2520	7	0.47	0.34	0.26	0.99
	Limari (LI)	11,650	37	1673	7	0.62	0.45	0.30	1.02
	Choapa (CHO)	7815	39	1701	7	0.46	0.33	0.29	1.02
	Aconcagua (AC)	7341	42	1847	7	0.37	0.37	0.24	1.20
	Average					0.47	0.36	0.26	1.06
NPS	Maipo (MP)	14,810	37	1664	8	0.49	0.38	0.25	1.28
	Mataquito (MAT)	6219	31	1106	7	0.20	0.22	0.18	1.13
	Rapel (RA)	14,041	33	1166	8	0.44	0.39	0.27	1.21
	Maule (MA)	14,788	17	432	8	0.57	0.34	0.24	1.17
	Itata (IT)	11,457	19	581	8	0.40	0.34	0.27	1.17
	Bio Bio (BB)	24,223	24	805	8	0.28	0.30	0.22	1.06
	Imperial (IM)	13,443	15	397	7	0.44	0.42	0.24	1.09
	Tolten (TO)	8100	22	555	7	0.39	0.35	0.20	1.15
	Valdivia${}^{*}$ (VA)	11470	23	489	7	0.49	0.39	0.28	0.96
	Bueno (BU)	13,897	19	422	8	0.46	0.40	0.30	1.12
	Palena (PAL)	11,584	41	865	8	0.49	0.20	0.24	1.00
	Aysen (AY)	12,781	36	834	7	0.62	0.35	0.29	1.01
	Average					0.44	0.34	0.25	1.11
AP	Baker${}^{*}$ (BA)	29,326	31	891	8	0.88	0.37	0.23	1.06
	Pascua (PA)	12,141	31	943	8	0.54	0.31	0.20	1.06
	Average					0.71	0.34	0.22	1.06

.

In order to explore the fractal dimension of the drainage networks at a more detailed scale, a set of 30 drainage sub-networks were extracted from Loa, Elqui, Valdivia and Baker basins (ordered in north–south direction). The same analysis was applied to each of these units obtaining the structure and morphometric properties of these networks. The features of these units can be observed in Appendix A. In this way, eight sub-networks were extracted from Loa network (referred as LO1 to LO8 in

Table 2. Morphometric parameters of large networks, where A is the drainage area, $i_m$ the mean slope, $H_m$ the mean elevation and $\Omega$ the maximum order of each network according to Strahler’s hierarchical ordering. The indices $F,C,E$ and the drainage density $\rho$, were included for reference.

Basin	Sub-Basin	A (km2)	${\mathit{i}}_{\mathit{m}}$ (%)	${\mathit{H}}_{\mathit{m}}$ (m)	$\mathbf{\Omega }$	F	C	E	$\mathit{\rho }$ (km${}^{-1}$)
Loa	LO1	120	6.17	974	4	0.24	0.23	0.56	1.21
	LO2	122	4.50	980	4	0.19	0.26	0.49	1.29
	LO3	7551	7.51	1985	8	0.36	0.25	0.68	2.21
	LO4	8017	8.70	3800	8	0.26	0.19	0.57	1.70
	LO5	311	4.89	1181	6	0.33	0.26	0.65	1.81
	LO6	616	8.30	1663	6	0.24	0.29	0.55	2.88
	LO7	469	5.57	1489	6	0.17	0.25	0.47	2.83
	LO8	3208	5.49	2170	7	0.37	0.21	0.69	2.51
	Average					0.27	0.24	0.58	2.05
Elqui	EQ1	1073	15.80	1027	6	0.32	0.23	0.64	1.05
	EQ2	737	18.81	1634	6	0.36	0.28	0.67	1.02
	EQ3	4086	25.48	3527	6	0.29	0.18	0.60	0.91
	EQ4	563	18.57	1052	5	0.26	0.20	0.58	0.97
	EQ5	261	22.16	1666	4	0.20	0.20	0.51	1.03
	EQ6	131	20.93	1368	4	0.21	0.29	0.52	1.23
	EQ7	51	21.42	1481	4	0.37	0.32	0.69	1.59
	EQ8	1515	27.15	3202	6	0.28	0.28	0.60	0.97
	Average					0.29	0.25	0.60	1.10
Valdivia	VA1	3367	8.31	215	6	0.15	0.18	0.43	1.09
	VA2	1486	13.76	630	6	0.19	0.22	0.50	1.10
	VA3	1386	18.81	997	6	0.37	0.16	0.68	0.99
	VA4	316	10.94	220	5	0.20	0.20	0.50	1.03
	VA5	615	9.29	200	5	0.28	0.23	0.60	1.01
	VA6	107	6.61	207	4	0.21	0.32	0.51	1.14
	VA7	960	11.76	447	6	0.29	0.23	0.61	0.11
	Average					0.24	0.22	0.55	0.92
Baker	BA1	1896	20.33	874	6	0.14	0.19	0.43	0.96
	BA2	3197	17.66	934	6	0.15	0.16	0.44	1.06
	BA3	394	27.62	1282	5	0.22	0.21	0.53	0.87
	BA4	306	14.68	898	5	0.09	0.28	0.34	1.00
	BA5	4785	8.32	870	7	0.40	0.31	0.72	1.16
	BA6	386	19.61	1061	4	0.30	0.25	0.62	0.91
	BA7	1499	18.82	1126	6	0.19	0.18	0.49	1.12
	Average					0.22	0.23	0.51	1.01

and

Table 3. Fractal dimension estimated from Fractalyse ($D_F$) for Loa, Elqui, Valdivia and Baker networks (marked in red). The networks were sorted by tectonic-segment influence. Sub-networks were also included. The mean error of each data, ${\Delta }_F$, was also included for reference.

Plate	Basin	Sub-Basin	${\mathit{D}}_{\mathit{F}}$	Plate	Basin	Sub-Basin	${\mathit{D}}_{\mathit{F}}$	Plate	Basin	Sub-Basin	${\mathit{D}}_{\mathit{F}}$
NPN	Loa		1.81 ± 0.02	NPS	Maipo		1.64 ± 0.05	AP	Baker		1.67 ± 0.05
		LO1	1.22 ± 0.05		Mataquito		1.50 ± 0.06			BA1	1.39 ± 0.06
		LO2	1.19 ± 0.04		Rapel		1.62 ± 0.06			BA2	1.45 ± 0.06
		LO3	1.66 ± 0.04		Maule		1.63 ± 0.06			BA3	1.21 ± 0.04
		LO4	1.62 ± 0.04		Itata		1.62 ± 0.07			BA4	1.24 ± 0.04
		LO5	1.33 ± 0.05		Bio Bio		1.65 ± 0.05			BA5	1.55 ± 0.06
		LO6	1.49 ± 0.06		Imperial		1.59 ± 0.06			BA6	1.27 ± 0.05
		LO7	1.46 ± 0.04		Tolten		1.57 ± 0.06			BA7	1.37 ± 0.05
		LO8	1.60 ± 0.04		Valdivia		1.56 ± 0.07		Pascua		1.62 ± 0.06
	Quebrada Caracoles		1.75 ± 0.03			VA1	1.44 ± 0.06
	Salado		1.66 ± 0.05			VA2	1.49 ± 0.04
NPFS	Copiapo		1.64 ± 0.05			VA3	1.37 ± 0.06
	Huasco		1.57 ± 0.06			VA4	1.24 ± 0.05
	Elqui		1.54 ± 0.06			VA5	1.29 ± 0.05
		EQ1	1.36 ± 0.05			VA6	1.18 ± 0.04
		EQ2	1.32 ± 0.05			VA7	1.36 ± 0.06
		EQ3	1.45 ± 0.05		Bueno		1.61 ± 0.07
		EQ4	1.31 ± 0.06		Palena		1.57 ± 0.06
		EQ5	1.23 ± 0.04		Aysen		1.58 ± 0.06
		EQ6	1.20 ± 0.03
		EQ7	1.22 ± 0.04
		EQ8	1.36 ± 0.06
	Limari		1.59 ± 0.07
	Choapa		1.55 ± 0.06
	Aconcagua		1.67 ± 0.04

), seven sub-basins were extracted from Elqui network (referred as EQ1 to EQ8 in Table 2 and Table 3); from Valdivia network seven more (referred as VA1 to VA7 in Table 2 and Table 3) and finally, seven sub-basins from Baker network (referred as BA1 to BA7 in Table 2 and Table 3). To the best of our knowledge, only Dorsaz et al. [47] have reported a detailed morphometric study for some Chilean networks, located between Pampa Colorada and Pampa Tamarugal. Our study extends such analysis by considering several units distributed across the continental territory and influenced by the tectonic segments acting along the Pacific coast.

3. Results

3.1. Geomorphological Description of Chilean Territory

According to detailed physical, geographical and geological description of Chilean territory provided by [48,49,50] and Moreno [51], it is possible to affirm that Chilean landscape is the result of endogenous and exogenous processes, acting simultaneously across the territory. Endogenous process generates from the subduction phenomenon of Nazca and Antarctic plates beneath the South American plate. This process was responsible for the morphostructural conformation of the western relief of extreme southern regions in South America. On the other side, exogenous processes arise due to subtropical and polar winds influence defining both, thermal and pluviometric regimes of the country. This process also defines the different climates developing in the country, contributing to modeling the morphostructural bands observed in Figure 1b. Three of such bands play a key role on drainage patterns later described: the Cordillera de los Andes, the Central Depression and the Cordillera de la Costa.

In this context, from 18${}^{\circ }$ S–27${}^{\circ }$ S morphostructural bands arrange in parallel. Andes mountains show very high altitudes of volcanic origin. The lithological types correspond to extrusive rocks, accumulating in regions of high elevation. In the zone known as Altiplano, also located in Andes mountains, the influence of subtropical Amazonian winds generate summer torrential precipitations. In this region, rivers transport sediment loads and water towards the central depression. Coastal mountains band show elevations lower than Andes relief. Such band separates from Central Depression by Atacama Fault Zone, extending around 1000 km in the north–south direction. Close the coastal line, coastal band appears in the form of cliffs showing in average altitudes of around 500 m. Except by Altiplano’s induced rainfall, the rest of this region shows hyper-aridity patterns. In such scenario, Atacama desert located in the Central Depression, the driest desert on earth, is crossed by some rivers showing very occasional runoffs. Towards the coast, fogs develop due to increasing evaporation because of the higher temperatures and Pacific Ocean influence. Nevertheless, such fogs do not penetrate beyond the coastal band. In this region, xerophyte-like vegetation arise, heterogeneously distributed across the region and adapted to low atmospheric humidity conditions. Low grasses, cacti and some shrubby species can be found. Soils layers are thin and very poor in organic matter content.

From 27${}^{\circ }$ S to 33${}^{\circ }$ S, the territory is dominated by mountainous landscapes in the east–west direction. Due to penetration direction of Nazca plate, Central Depression has almost disappeared and volcanic activity is almost nonexistent. Relief rises in the west–east direction and minor fault systems can be recognized. These characteristics avoid a clear distinction between Coastal and Andes mountains. Morphostructural organization looks fragmented, mainly by deep fluvial valleys showing an east–west orientation. This region is governed by less arid climatic conditions leading to a clear transition between hyper-arid climate in the north and tempered-zones in the south. Vegetation and soils layers are well developed across the valleys due to significant supply of nutrients and the presence of long runoffs, both of them essential for agriculture. Shape of transverse valleys help coastal haze to penetrate into continental territory, providing moisture conditions for vegetation development. Such vegetation does not form, however, a homogeneous cover.

From 33${}^{\circ }$ S to 42${}^{\circ }$ S, morphostructural bands arise again. Andes mountains band shows active volcanism and Central Depression is bounded by Andean landscapes and Coastal mountains. However, such reliefs show lower elevations. Here, temperate climate becomes more humid and is seasonally influenced by polar fronts. Therefore a seasonal contrast can be well distinguished across this region. Summers are dry, hot and long. Meanwhile, winters are cold and wet. Under these conditions, soil layers are deep allowing the development of important agricultural activities. Such activities are essentially conducted into the Central Depression supported by abundant runoff. Vegetation is essentially of mesophyte type, dominated by sclerophyllous scrub Finally, from 42${}^{\circ }$ S to extreme south, landscape is fragmented due to Quaternary glaciation process practically covering the whole territory in the east–west direction. Central depression disappears submerging into the ocean, as well as, coastal mountains can be observed until the latitude 46${}^{\circ }$ S. In this last region, Taitao triple junction arises. In this zone Nazca, Antarctic and South American plates join. Towards the south, seismic and volcanic activity is less frequent, a characteristic closely related to moderate subduction speeds. Andes mountains are mostly formed by granitic rocks (batholiths), that can reach the ocean and formed by many small islands, with block tectonics and significant structural faults. Here, the Liquiñe-Ofqui Fault Zone together with erosive glaciation effects can be observed. Climate conditions are extremely cold and rainy during the year. Vegetation has enough water availability to develop during the year and is dominated by forests covering the region.

3.2. About the Drainage Patterns of Chilean River Networks

Figure 2 shows the river drainage patterns obtained for large networks. In general, the geometry of each network is the combined result of slope effects and landscape organization, all of them described in Section 3.1. Concerning this last point, dendritic patterns become the dominant feature in Chilean river networks as observed in Figure 2. Drainage patterns also show a clear east–west runoff direction and they are generated over relatively homogeneous lithological substrates, showing similar resistance to hydrodynamic erosion with tributaries connecting at acute angles (<90${}^{\circ }$). At a more detailed scale, some units present grid or rectangular drainage patterns a sign of structural control on streams diffusion. The largest river networks develops over the three aforementioned morphostructural bands. Therefore, characterizing the topology of these drainage networks just by looking at their patterns and drainage density is not enough, justifying the use of fractal dimension to take into account such variability.

According to Table 1, Loa, Caracoles and Salado networks in the Nazca North plate observed in Figure 2 show high drainage densities, that is, $\rho =1.86,2.15,1.55$ km${}^{-1}$, respectively. For these units the form factor is quite variable too ranging from 0.28 to 1.88. Basin circularity and elongation indices seems to be, however, uniformly distributed. Following the description given in Section 2.2, these last parameters suggest the basins are tectonically activated dominated by permeable soils. The average slopes of these units also falls in the range $9\%\le i_m\le 20\%$, which are particularly high. The northern networks develop mainly over the Central Depression in a hyper-arid climate context and under the tectonic influence of the northern zone of the Nazca plate. Their drainage patterns are essentially dendritic although, in Salado network the influence of the Central Depression begins to disappear. Here, part of the basin develops on mountainous sectors and rectangular patterns are also locally observed. At sub-basin level, networks LO1 to LO8 in Table 2 show instead lower values of the form factor, a more homogeneous distribution of the circularity index (C < 0.29) and the elongation index ranging from 0.49 to 0.69, higher than measures reported for Loa network in Table 1. This is a sign of elongated shapes dominated by steep gradients, high elevations and slight tectonic activity, a natural feature of the landscape in those regions of the country.

Respect to their drainage patterns and density, networks located between Copiapo and Aconcagua rivers look very similar. Drainage density for these units ranges from 0.98 to 1.20 km${}^{-1}$, lower than values observed for Nazca North plate. Many factors could be influencing this parameter, e.g., increase of vegetation and rainfall intensity. These networks develop mainly over mountainous landscapes where Central Depression has started to disappear, falling into the influence of Nazca Plate flat subduction (Nazca Flat Slab). According to Table 1, the form factor here ranges from 0.37 to 0.62, basin circularity $C\le$ 0.45 and $E\le$ 0.30. Values of E correspond to elongated basins with active tectonic influence. According the description provided by Moreno [51], such influence is a common feature for large networks explored in this study. These last systems organize around a main stream and the dominant drainage pattern is once again dendritic, especially in those regions of Central Depression. However, there is a clear trend to show lattice patterns in mountains regions of Cordillera de los Andes and Cordillera de la Costa. There, smaller streams connect to larger ones at almost 90${}^{\circ }$. These characteristics traduce in quite homogeneous drainage densities fluctuating in the range 0.98 km${}^{-1}$$\le \rho \le$ 1.20 km${}^{-1}$ and 0.37 $\le F\le$ 0.62. Average slopes for these units are particularly high with values falling in the range 33% $\le i_m\le$ 46%. Same data analyzed at sub-networks level (EQ1 to EQ8 in Table 2), show in average that $E\approx 0.60$ (slightly higher than values reported for Loa sub-networks), $C\le$ 0.32 and F = 0.29. These last values are quite different from those obtained for Elqui network as observed in Table 1.

Basins between Maipo and Aysen rivers, locate to the south of Nazca Plate (referred as NPS in Table 1). They develop under the influence of mediterranean and tempered climates. In this region, drainage pattern and density show a clear influence of the aforementioned morphostructural units. This influence traduces in north–south strips arrangement. In almost all the cases, a dendritic drainage pattern is observed again over Central Depression and mountain regions. Lattice and rectangular patterns can also be observed with tributaries bifurcating at almost 90${}^{\circ }$. It is interesting that Itata network shows a parallel pattern developing very well across the Central Depression. This case corresponds to large fluvio-alluvial fans associated with Laja river sub-basin, formed from successive dam breaks of Quaternary volcanic debris. Here, F ranges from 0.20 to 0.62 and $\rho$ from 0.96 km${}^{-1}$ observed in Valdivia network (close to the southern limit of this plate), to 1.28 km${}^{-1}$ obtained for Maipo network (in the extreme north of this plate). The parameters C and E remains below 0.42 and 0.30, respectively. However, sub-networks VA1 to VA7 shows in average that F = 0.24, C = 0.22 and E = 0.55. Even though these sub-networks locate in a different tectonic environment, their values are not significantly far from those observed for LO1 to LO8 and EQ1 to EQ8, leading to almost homogeneous shapes at that scale of observation.

Finally, Palena and Aysen basins located in Austral and Patagonian regions, fall under the tectonic influence of Antarctic plate developing over Cordillera de Los Andes. Common drainage patterns here are parallel and trellis. Acute contact angles between tributary streams are practically nonexistent. For this reason, the structural control of these networks is essentially linked to crust intense faulting and aggressive erosive processes dating from Quaternary glaciation period. The temperate climate of this region provides abundant rainfall conditions throughout the year. These characteristics strongly influence rivers runoff regime, many of them showing torrential high flow rates. For these networks the average slope is $i_m$≈ 31% (very high) and drainage densities are particularly high, falling in the range 1.87 km${}^{-1}$$\le \rho \le$ 2.49 km${}^{-1}$. At sub-basin level, the morphological parameters of each network (BA1 to BA7) are in average F = 0.22, C = 0.23 and E = 0.51 as shown in Table 2. Once again, at a sub-network level no dramatic differences can be observed between these averaged values respect to those reported above. However, these morphometric indices do show a significant departure from those measured at large networks level.

3.3. Distribution of Horton Ratios

Figure 3a–c show the frequency distributions of the average Horton ratios for all the networks (23 large networks and 30 sub-networks). A first look shows that peaks of frequency are reached at different ratios, that is, $R_A=5.00,R_B=4.25,R_L=2.75$ for large networks and $R_A=5.75,R_B=5.25,R_L=3.25$ for sub-networks. These values are close to some measurements reported in the literature [34,35,52]. Notice also that distributions observed in Figure 3a–c are not qualitatively different. However, distribution calculated for all networks shows a significant departure from the rest of the curves (Figure 3d). This last observation can be also observed in Figure 4, showing the Cumulative Normalized Frequency Distribution (CDF) of Horton ratios, while large units show a quick saturation, sub-networks evidences a more extended distribution, then a size effect can be observed in our measurements. An arbitrary value CDF = 80% was included to emphasise this effect (horizontal dashed line in Figure 4a–c). Notice that distributions related to $R_A$ and $R_B$ are not very different, particularly those related to large networks. However, distribution related to $R_L$ shows a significant departure respect to the rest of Horton parameters. This behaviour holds when looking at the joint distribution (Figure 4c). In particular, two different behaviours of each distribution can be observed around when CDF ≈ 80%. When CDF > 80% the shape of these distributions is essentially controlled by measurements conducted at sub-network level.

Figure 5 shows the distribution of $R_A,R_B,R_L$ versus $log\left(A\right)$ (the logarithm of network area) for all units. Measurements were splitted-up into two datasets: sub-networks, in the range $1.7\le log\left(A\right)\le 3.8$ and large networks, in the range $log\left(A\right)>3.8$. The error bands, ${\Delta }_B,{\Delta }_L,{\Delta }_A$, were also included, and estimated as the mean quadratic error of the slope of the curves $log\left(N_{\omega }\right),log\left({\overline{L}}_{\omega }\right),log\left({\overline{A}}_{\omega }\right)$ versus $\omega$, respectively. These distributions were also provided by GRASS-GIS. In our study, these bands varied in the ranges 0.36 $\le {\Delta }_B\le$ 2.07, 0.27 $\le {\Delta }_L\le$ 1.75 and 0.28 $\le {\Delta }_A\le$ 2.89 for sub-networks and, 0.23 $\le {\Delta }_B\le$ 0.65, 0.13 $\le {\Delta }_L\le$ 0.35 and 0.15 $\le {\Delta }_A\le$ 0.52 for large networks. Arbitrary scattering bands of $\pm 50\%$ around each fit were drawn for emphasizing measurement dispersion. Surprisingly, each distribution shows a decreasing behaviour that can be reasonably fitted by a linear function of slope $\eta <0$, that is, $R_q=\eta log\left(A\right)+\tilde{c}$ as observed in Figure 5, where $\eta ,\tilde{c}$ are free fitting parameters. The slope is quite similar for the curves $R_L$ and $R_B$ versus $log\left(A\right)$ ($\eta =-$0.78 and −0.73, respectively). However, we obtain $\eta =-$1.40 when comparing $R_A$ and $log\left(A\right)$. This last result indicates that ratio $R_A$ decays faster than the rest of Horton’s parameters.

3.4. Fractal Dimension of Chilean Networks

In this section, the interplay between the fractal dimension and the morphometric properties of networks is explored. Table 3 shows the fractal dimension of each network obtained from Fractalyse, ordered by the influence of tectonic segments. The error band of each data, ${\Delta }_F$, was included for reference. If we look at large networks dataset, we observe that those units located at Nazca North plate-segment (referred as NPN in Table 3) show the highest values of the entire record. In particular, the fractal dimension of Loa network is $D_F\approx 1.81\pm 0.02$, the highest dimension of the present study. However, when moving towards the Nazca Flat Slab plate (referred as NPFS), the fractal dimension of the units decreases to values in the range 1.50 ± 0.06 to 1.67 ± 0.04. By contrast, units influenced by Nazca South plate (referred as NPS) show values distributed ranging from 1.56 ± 0.07 to 1.61 ± 0.07, similarly to those networks located in front of NPFS plate segment. However, the fractal dimension of networks of Austral-Patagonian regions (referred as AP) increases again, leading to values in the range 1.62 ± 0.06 to 1.67 ± 0.05 corresponding to Pascua and Baker networks, respectively. Thus, significant differences arise on the empirical fractal dimension depending on the tectonic segment influencing the landscape. These differences suggest that the parameter $D_F$ could be controlled by a competition between tectonic and erosive-related processes, following the interpretation proposed by Phillips [31] and Donadio et al. [33].

Figure 6a compares the parameter $D_F$ with the logarithm of the mean slope log$\left(i_m\right)$ of the networks. Sub-network dataset was also included for comparison. Not clear correlations can be deduced between these points, however, most of large networks measurements fall into a moderated to high-slope regions. Loa’s sub-networks (LO1 to LO8) show, by contrast, lower mean slopes than the basin containing them (indicated as LO in Figure 6a,b). A similar observation is obtained when analysing sub-networks from Elqui, Valdivia and Baker. This suggests that some morphometric properties analysed at a more detailed scale of observation (also called “mesoscale”) show a significant departure from the properties observed at larger scales. This is a first sign of a behaviour not consistent with the hypothesis of self-similarity. A large subset of measurements shown in Figure 6a falls into a moderated slope region, constituting a transition-like region between low and high-sloped data. Figure 6b compares $D_F$ with the drainage density $\rho$ of each network. Most of points concentrate into the band defined by 0.8 km${}^{-1}$$\le \rho \le$ 1.24 km${}^{-1}$, except by sub-network VA7. In this region, $D_F$ ranges from 1.20 ± 0.03 to 1.67 ± 0.04 and once again, no reasonable function seems to fit this data. When $\rho >1.24$, the fractal dimension shows a slight increase for increasing values of the drainage density leading to the value $D_F=$1.81 ± 0.02 obtained for Loa network. Notice that into this last range of density, most of data correspond to networks tectonically influenced by Nazca North plate-segment.

Figure 7 compares $D_F$ with $log\left(A\right)$. Three well ordered datasets arise from this comparison. Each dataset can be characterized by a linear fit of slope s, that is, $D_F=s·log\left(A\right)+\xi$, with $\xi$ a free fitting parameter. In this study, we have obtained $s\approx 0.33$ for most of large basins, $s\approx 0.25$ for sub-networks EQ7, LO6 and LO7 and a central fit of arbitrary slope $s\approx 0.30$, for the rest of units. This organisation pattern is quite unexpected considering the disorganised character of the Chilean territory with respect to their latitudinal development, the current and past tectonic and the geomorphological processes and the different factors conditioning the climatic and hydrological characteristics of the basins already discussed in Section 3. The fits proposed in Figure 7 suggest then that $D_F$ depends upon network’s area, a feature that was reported in urban environments [53].

Interested on exploring the differences between the empirical results provided by the box-counting method and the analytical models proposed in Section 1, in Figure 8 we compare $D_F$ with the models of La Barbera and Rosso [2,24], Tarboton et al. [21] and Liu [20]. Error bands of $D_1,\cdots ,D_4$ and $D_F$ were included for reference. Such bands corresponds the range of variation of ratios $\frac{log(R_B\pm {\Delta }_B)}{log(R_L\pm {\Delta }_L)}$ and $\frac{log(R_B\pm {\Delta }_B)}{log(R_A\pm {\Delta }_A)}$. Although data dispersion is non negligible, equation $D_1$ proposed by [2] shows reasonable correlation with Fractalyse results, for both kind of networks. By contrast, equation $D_4$ established by [20] overestimates this dimension leading to the highest values of our record. For some basins $D_4$ is higher than 2, judged unrealistic given the lower drainage densities of our networks. Then this last model might be used with caution. Models proposed by [24] and [21], $D_2$ and $D_3$, respectively, provide results strictly comprised between the results obtained from $D_1$ and $D_4$. Two linear fits were calculated to emphasize this observation, that is, $D_1=0.94D_F$ and $D_4=1.17D_F$ in Figure 8a and $D_1=0.98D_F$ and $D_4=1.31D_F$ in Figure 8b. The Pearson correlation coefficient $r_k$ was also estimated between $D_k$ and $D_F$, for $k=1,\cdots ,4$ (following the previous notation). The following values were obtained $r_1=r_2=r_3=0.81$ for large networks and, $r_1=r_2=r_3=0.83$ for sub-networks, showing high correlation between data. In the rest of cases, a poor performance of such index was obtained ($r_4=$0.31 and 0.47, respectively). Then a reasonable agreement between $D_1$ and empirical measurements provided by the box-counting method, can be established.

To finish this section we explore how far both methods, Fractalyse and Horton’s law, could be. For this purpose the parameters $D_F$ and $D_k$ were compared, for $k=1,\cdots ,4$. The influence of drainage area A on this possible connection was also explored. First, notice that ratios $log\left(R_B\right)/log\left(R_L\right)$ and $log\left(R_B\right)/log\left(R_A\right)$ are involved on the dimensions $D_k$. Then, for convenience, the auxiliary parameter $\zeta =log\left(R_p\right)/log\left(R_q\right)$ was defined, where $p=B$ and $q=L$ or $q=A$. The relationship between $\zeta$ and A is shown in Figure 9a. Referential averages were included for each distribution, that is, $\overline{\zeta }=1.41\pm 0.14$ for $q=L$ and $\overline{\zeta }=0.92\pm 0.04$ for $q=A$. When $q=L$ (points in red) measurements show a significant departure respect to its average, particularly in the range $A\ge$ 6219 km${}^2$. Such departure is, however, negligible when $q=A$ (points in gray), leading to measurements almost independent of drainage area. Thus, a relationship between $\zeta$ and A arises depending on the model used for computing $\zeta$. Figure 9b compares the ratio $D_F/\zeta$ and $log\left(A\right)$, for $q=L$ and $q=A$. Notice that points in gray show significant dispersion in the whole range. When $q=L$, the ratio $D_F/\zeta$ is slightly larger than 1, showing a linear increase with $log\left(A\right)$ for $A<$ 6219 km${}^2$. Such increase is then, followed by a saturation-like effect for very large drainage areas ($log\left(A\right)\ge 3.79$). In the second case ($q=A$), the ratio $D_F/\zeta$ is almost constant. Both curves suggest that Fractalyse and Horton’s law approaches lead to close results when dealing with large networks, except for a constant value. However, strong differences arise when analysing small-networks, where an effect of drainage area can be observed.

Then, from the distributions shown in Figure 9 data related to $q=L$ can be roughly described by the fitting function $D_F/\zeta ={\kappa }_{1}log\left(A\right)+{\kappa }_2$, which is valid in the range A < 6219 km${}^2$. Meanwhile, we obtain $D_F/\zeta ={\kappa }_3$ in the range $A\ge$6219 km${}^2$. The fitting coefficients in both regimes are $({\kappa }_1,{\kappa }_2,{\kappa }_3)=(0.89\pm 0.05,0.04\pm 0.02,1.06\pm 0.01)$. By contrast, when $q=A$, we obtain $D_F/\zeta \approx {\kappa }_4$ for the whole range of A with ${\kappa }_4=1.13\pm 0.01$. Therefore, a significant disagreement between both approaches arise when considering $R_L$ or $R_A$ on these computations. Such differences emphasize the drainage area effect on determining the fractal dimension. However, according to Figure 9 this effect is important only into the range $A\ge$ 6219 km${}^2$. Beyond this point, both approaches conduct to dimensions almost independent on A separated for a constant. Third, the saturation effect observed in Figure 9b re-affirms the idea that the fractal dimension of a river network cannot growth indefinitely. An asymptotic limit must arise for these distributions in agreement with the full-filling space concept proposed by Phillips [31], which imposes that fractal dimensions cannot exceed the topological dimension of the embedding space (=2) for two-dimensional objects.

In conclusion, the distributions observed in Figure 9b suggest for this study that a relationship can be established between the fractal dimension obtained from Horton’s law and Fractalyse. Such relationship is of the type $D_F=\psi (A,d)D$, where d is the meandering fractal dimension of mainstreams introduced by [22], D the fractal dimension of the whole network based on Horton’s law (=$D_1,D_2,D_3$ or $D_4$ in this article) and $\psi (A,d)$ a function whose structure must be established according to regional characteristics of basins. In this paper, we show that at sub-network level the function $\psi (A,d)$ describes a logarithmic variation with A, although for large networks we obtain $\psi \approx c$, where c is a given constant.

4. Discussion

4.1. Chilean River Networks: Self-Similar or Self-Affine Structures?

Considering all the results already shown, in this section we discuss to which extent Chilean networks can be treated as self-affine or self-similar objects, by following the definition given in Section 1. To this concern, notice that the values of $D_F$ already evidence a clear contrast acroos the different regions of the country, particularly for extreme regions (e.g., Loa vs. Baker regions). Meanwhile, river networks influenced by Nazca Flat Slab and Nazca South segments show more uniform values particularly, those units developed over the alluvial fans of the morpho-structural band called Central Depression (shown in Figure 1b), all of them influenced by Mediterranean climate. These latitudinal differences emphasize, once again, that structural controls might exist on these systems which seem to be behind the morphological anisotropy of Chilean landscapes. Such anisotropy constrains the diffusion of streams in agreement with Phillips [31]. This effect is present in most of our networks. Although some cautions have been raised about the extended use of the box-counting method to estimate the fractal dimension of real networks [40,54], in this study we show that this method provides reasonable values describing the morphological inhomogeneity of Chilean river systems.

Morphological, geological, tectonic and climatic processes acting across the continental territory, all of them described in detail in Section 3.1 and Section 3.2. Such processes contribute to the anisotropy of landscapes, as well. These characteristics give rise a clear east–west runoff pattern, with decreasing mean slopes of basins in the same direction, facilitating the existence of preferential directions on stream diffusion. Vegetation, lithology and urban settlements can be considered as second-degree factors, contributing to limiting this sculpting process and influencing the final drainage pattern described in Section 3.2. These agents, acting together with climatic and hydrological factors, play a key role on shaping the final geometry of drainage networks. These effects become notorious when moving towards the south of the country. These patterns of spatial variability define a self-affine character for Chilean river networks, implicit in Horton’s ratios as observed in Figure 10. In this study, very well ordered distributions arise for both kind of networks, all of them fitted by a power-law function:

$$R_p={R_q}^m,wherep\ne q,$$

(7)

However, exponents m are very different in both cases. At sub-network level, we obtain $m=1.26\pm 0.02$ when $p=B,q=L$ (see curve $C1$) and, $m=1.13\pm 0.01$ when $p=A,q=B$ (curve $C2$). However, at large network level, we obtain $m=1.49\pm 0.01$ when $p=B,q=L$ (curve $C3$) and, $m=0.98\pm 0.01$ when $p=A,q=B$ (curve $C4$). In consequence, the distributions $R_p=f\left(R_q\right)$ do not collapse into the same curve. This observation holds even when small-networks are extracted from the same basin. This last result suggests that streams networks morphology might be different at both observation scales. However, anisotropy or structural control effects are not possible to infer just by analysing Horton’s law. These ratios fail for such a purpose, because they do not take into account the shape of individual streams of a given network [40]. Maybe, this is the main difference respect to box-counting method, justifying the need of using both approaches to describe these kind of systems.

To emphasize this last observation, the fractal dimension of sub-networks was averaged and later compared to the dimension of the basin from which they were extracted. Among the different ways for computing such average (e.g., arithmetic mean), a weighted area-average better takes into account the spatial representativeness of fractal dimension calculated for each sub-network. Then, assuming a given network formed by n sub-networks, this weighted average $\overline{D}$ can be calculated as follows:

$$\overline{D}=\frac{D_1{A}_1+\cdots +D_n{A}_n}{A_1+\cdots +A_n},$$

(8)

where $A_i$ is the area of each sub-network of fractal dimension $D_i$, with $i=1,\cdots ,n$. The parameters $D_i$ were estimated from the models proposed by [2,20,21,24]. Loa, Elqui, Valdivia and Baker networks were chosen for this comparison. The results can be observed in

Table 4. Averaged fractal dimension $\overline{D}$ obtained from Equation (8) for Loa, Elqui, Valdivia and Baker, calculated according to models $D_1,\dots ,D_4$. Fractalyse values ($D_F$) were also included for comparison: ${}^{(*)}$ the averaged value of sub-basins and ${}^{(**)}$ the dimension of the respective large basin extracted from Table 3.

Weighted Average	When ${\mathit{D}}_{\mathit{i}}$ Is Calculated from the Model of:					From Table 3:
$\overline{\mathit{D}}$	[2]	[21]	[24]	[20]	${{\mathit{D}}_{\mathit{F}}}^{(*)}$	${{\mathit{D}}_{\mathit{F}}}^{(**)}$
Loa	1.58 ± 0.01	1.80 ± 0.11	1.84 ± 0.11	1.85 ± 0.04	1.61 ± 0.04	1.81 ± 0.02
Elqui	1.35 ± 0.04	1.54 ± 0.04	1.57 ± 0.04	1.84 ± 0.07	1.39 ± 0.05	1.54 ± 0.06
Vadivia	1.35 ± 0.07	1.54 ± 0.08	1.57 ± 0.08	1.81 ± 0.04	1.41 ± 0.05	1.56 ± 0.07
Baker	1.36 ± 0.05	1.55 ± 0.06	1.58 ± 0.06	1.86 ± 0.04	1.45 ± 0.05	1.67 ± 0.05

. Table 4 evidences serious discrepancies respect to the model used for computing these averages. For example, averaged Fractalyse values at sub-network level are far from data presented in Table 3. A similar conclusion can be obtained for the rest of networks and adopted models. Then, there is no way of recovering the fractal dimension of a large basin by simply averaging the fractal dimension of its portions. Thus, self-similar patterns of diffusion are difficult to observe in Chilean fluvial networks by far, as suggested by Kim [15,32]. A quantitative determination of the self-affinity degree of our networks could help to discriminating these effects, for example by computing the Hurst parameter [23].

To finish this section, we point that the role played by the mean slope and the drainage density of the basin remains unclear as shown in Figure 6. Data collected for small-networks organize in ranges different from those observed for larger networks as observed in Loa, Elqui, Valdivia and Aysen basins. These changes violate the scale-invariance principle expected for self-similar objects, an observation that has been supported by Kim [32]. Unexpectedly, when $D_F$ is compared with the drainage area clear groupings can be observed, described by the scaling law $D_F\propto s·log\left(A\right)$, where s is an area-dependent exponent (see Figure 7). A similar result was reported by Shen [53] from measurements conducted in urban environments. In consequence, tectonics and hydrological processes influence the final drainage patterns acting in different ways for large and sub-networks. Both processes are the key to understand how Chilean drainage patterns evolves in time, as suggested by Turcotte [55].

4.2. A Relationship between $\rho ,i_m,D_f$ and Tectonics

From the previous results, we can readily analyse the relation between the mean slope, the drainage density and the fractal dimension $D_F$ of each network. Such relationship can be observed in Figure 11. Tectonics effects were also qualitatively indicated in the same diagram, where we can remark two things. First, the subduction influence of Nazca and Antarctic plates shows a clear west–east direction, also known as the collision margin with the South American plate. Such phenomenon has generated the main morphostructural units organized in longitudinal strips of the relief explained in Section 3. Second, in general the development of the different climatic types is related with the global atmospheric circulation patterns, inducing latitudinal climatic zones in the north–south direction. For this reason, arid climates in the north progressively becomes more humid towards southern regions. From these considerations, it is reasonable to think about isolated or combined effects in the development of the geometrical patterns of drainage networks.

In this context, most of river networks in Chile have been the combined result of tectonics and erosive processes, giving rise to very coarse drainage textures, showing moderate to high slopes in large networks. For these last units, the mean values of the fractal dimension $D_F$ ranges from 1.66 to 1.81, with networks falling under the influence of Nazca Plate-North (NPN). These values show a quick decrease to the south of the country in agreement with the change on tectonics and the power of erosive process. The highest values of the fractal dimension of our record (Loa, Quebrada Caracoles and Salado rivers) are also observed in regions dominated by moderate to low slopes, with coarse drainage textures ($\rho >1.24$). Most of networks located in Nazca Plate-Fat Slab (NPFS) and Nazca Plate-South (NPS) show a quasi-uniformity on the values of the fractal dimension. Coincidentally, most of these units show very coarse drainage texture with a wide range of variability of its mean slope. This suggests that the fractal dimension is not lonely conditioned by the slope of the relief, but also another effects begin to control the drainage process (e.g., sediment transport in rivers, landslides, erosion induced by rainfall, vegetation effects). This scenario changes again for networks located at the Antarctic Plate (AP), where the fractal dimension increases again, ranging from 1.62 to 1.57 for large networks and from 1.21 to 1.55 for sub.networks. This increase on $D_F$ could be explained a priori by the change on the erosion mechanisms involved on the drainage process (e.g., glacier erosion and deposition). We believe that this level of organisation is not coincidental and reveals a non-aleatory organization between tectonics, climate regimes and fractal morphology, questioning the apparent disordered character of the territory mentioned by Subercaseux [7].

5. Conclusions

In this paper, we have measured for the first time the fractal dimension of Chilean river networks, analysed at two different size scales. Chilean networks are complex geometrical objects that show variable patterns of ramification, inserted into landscapes of different morphological characteristics. In such a context, the fractal dimension estimated from the box-counting method, computed through the software Fractalyse, helps to discriminating the drainage patterns observed across the territory, from northern to southern regions. Such dimension also reveals a significant contrast between large and sub-networks, even when these last ones are extracted from the same morphological unit. The box-counting method captures not only the overall properties about the structure of the network, but also the individual properties of the streams. This establish by one hand, a significant advantage respect to the analytical methods based on Horton’s law and, on the other side, it allows to discriminate if a structural control is present in these systems. Then some cautions need to be considered when using analytical models based for this purpose as underlined by [17,26,56].

Despite these differences, both approaches still provide rich quantitative information about the geometric characteristics of river systems, compared to the information extrapolated from classical morphometric indices such as the form factor or the circularity index introduced before. In this way, fractal dimension even calculated as a single number for a drainage network, go farther establishing some connections with the sculpting process of the landscape as demonstrated by [15,23,32,33,34,35,52,57,58,59] and recently by Feng [36]. Certainly, accurate estimations of the fractal dimension will improve our observations which requires more sophisticated methodologies. To this concern, a multifractal analysis arises as a powerful tool that could capture the local characteristics of ramification patterns and the main spatial properties involved in Chilean river networks, across different scales.