Singularity cities

Abstract

We propose an upgraded gravitational model which provides population counts beyond the binary (urban/non-urban) city simulations. Numerically studying the model output, we find that the radial population density gradients follow power-laws where the exponent is related to the preset gravity exponent γ. Similarly, the urban fraction decays exponentially, again determined by γ. The population density gradient can be related to radial fractality and it turns out that the typical exponents imply that cities are basically zero-dimensional. Increasing the gravity exponent leads to extreme compactness and the loss of radial symmetry. We study the shape of the major central cluster by means of another three fractal dimensions and find that overall its fractality is dominated by the size and the influence of γ is minor. The fundamental allometry, between population and area of the major central cluster, is related to the gravity exponent but restricted to the case of higher densities in large cities. We argue that cities are shaped by power-law proximity. We complement the numerical analysis by economics arguments employing travel costs as well as housing rent determined by supply and demand. Our work contributes to the understanding of gravitational effects, radial gradients, and urban morphology. The model allows to generate and investigate city structures under laboratory conditions.

Keywords

Gravity models population density urban fraction fractal geometry

Introduction

The large number of processes working in cities make them complex objects extending over a range of spatio-temporal scales (Barthelemy, 2016; White et al., 2015). As pointed out by Batty (2013), a city science that explains city growth, sprawl, etc. needs to be supported by theories about how people relate to each other. Despite ongoing digitalization and globalization, geographical proximity still matters (Morgan, 2004). The small distances within cities, as extreme agglomerations, attract urbanites and thereby enhance the proximity.

Certainly, ideas of geographical gravitation have a long tradition and can be traced back to the middle of the 20th century and beyond (Carrothers, 1956; Stewart, 1948; Zipf, 1946). In view of new empirical findings, we revisit and extend a probabilistic city model (Rybski et al., 2013) from two states (non-urban, urban) to population counts. Specifically, we validate it against recent findings of urban fraction and population density gradients (Lemoy and Caruso, 2020) as well as of building heights within cities (Schläpfer et al., 2015).

The model to a large extent reproduces the features described for real-world cities. The numerical simulations enable us to relate both works to each other as well as to other properties including four different measures of city fractality and the fundamental allometry, i.e. between population and area of cities. Interestingly, the population density gradient decaying with the radial distance to the power −2 as found in (Lemoy and Caruso, 2017) corresponds to a fractal dimension of 0, which supports the point character of cities, i.e. singularities in space. We complement the numerical analysis by economics arguments employing travel costs as well as housing rent determined by supply and demand.

Model

We consider a two dimensional square lattice of size N × N whose sites w_j with coordinates $j \in {(1 \dots N, 1 \dots N)}$ can be empty (0) or occupied by an integer number of ‘inhabitants’. The probability that an inhabitant is added to a site is

q_{j} = G \frac{\sum_{k \neq j} w_{k} d_{j, k}^{- γ}}{\sum_{k \neq j} d_{j, k}^{- γ}}

(1)

where

d_{j, k}

is the Euclidean distance between the sites j and k. The denominator compensates for border effects of the finite system. The exponent

γ > 0

is a free parameter that determines how strong the influence of occupied sites decays with the distance. The constant G (exogenously) determines the overall growth rate and is given by

G = g / (max \frac{\sum w d^{- γ}}{\sum d^{- γ}})

where the parameter g is used to tune the speed of growth (

0 < g \leq 1

We start with an empty grid (w_j = 0 for all j) and, without loss of generality, put one inhabitant on the single central site. In every iteration, a random number z is drawn (from a uniform distribution between 0 and 1) for each grid cell with coordinates j and if $z < q_{j}$ then w_j is incremented by 1. We consider w as population counts in each grid cell. The procedure is repeated and stopped before the major central cluster reaches any of the system boundaries. Figure 1 shows examples of the emerging structures.

Figure 1.

Examples of city structures generated by the model, equation (1). For panels (a) and (b) $γ = 2.25$ was used and for (b) and (c) it was $γ = 2.5$ . For better comparability, in (a) and (c) as well as (b) and (d) iterations with approximately the same fraction p of occupied cells were chosen. The color bar indicates the number of ‘inhabitants’ in the cells (in log-scale). The red line in (a)-(d) shows the boundary of the major central cluster. Panel (e) is a 3-dimensional illustration of the major cluster from (d), whereas the third dimension is in linear scale.

This version differs from the original model (Rybski et al., 2013) only by (i) the w_j which originally were 0 or 1 and (ii) the g which originally was fixed to g = 1, so that the maximum probability was 1 (see Rybski et al., 2013, for details).

For some analyses, we extract the major cluster by applying the City Clustering Algorithm (CCA) (Fluschnik et al., 2016; Hoshen and Kopelman, 1976; Kriewald et al., 2016; Rozenfeld et al., 2008, 2011) with l = 1, i.e. only connecting nearest neighbors. The area $A_{c}$ of the major cluster is given by the number of cells with w > 0 belonging to the cluster. Analogously, the total population $S_{c}$ of the cluster is defined by the sum over w it consists of. For each major cluster, we extract its envelope, i.e. those cells which have at least one empty (nearest) neighbor which is not part of a hole within the cluster. We denote the number of cells the envelope consists of as perimeter C, and the largest distance from the envelope to the central cell of the lattice as $R_{c}$ .

The cells of the grid can also be understood as plots for buildings and the w_j as the height of the buildings. Assuming each floor corresponds to one apartment and each apartment is home to one person, then w_j corresponds to the number of inhabitants. More apartments per floor or more persons per apartment only represent a factor. We assume homogeneity, i.e. living space per person is constant throughout the city.

Analysis

On a square lattice of size 1000 × 1000 we run 10 realizations for various γ-values. As with the same normalization constant g a larger γ requires more iterations to fill the lattice, we take different normalization constants g for different γ-values to balance between the need of enough iterations and the computational time. Specifically, we run simulations for $γ = 2.0, 2.05, \dots, 2.7, 2.75$ with $g = 0.02, 0.02, 0.1, 0.1, 0.1, 0.1, 0.2, 0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1, 1$ , respectively. All iterations where the major central cluster is smaller than 200 occupied cells are excluded during the post-processing. We end up with approximately 1700 useful iterations in total for $γ = 2.0$ . The number of iterations increases with γ.

Radial gradients

First we want to study the gradients generated by the model and compare them with empirical results (Guérois and Pumain, 2008; Lemoy and Caruso, 2020; Peiravian et al., 2014). Following Lemoy and Caruso (2020) we define concentric rings around the initial central cell and calculate within them the population density and urban fraction. We also apply the rescaling proposed in Lemoy and Caruso (2020).

Population density gradient

The density is given by $D (r) = \sum_{k} w_{k} / \sum_{k} 1$ , where k is the index representing all cells at a distance between r and $r + δ r$ from the center. In this study the width of the rings is $δ r = 1$ . We only take rings up to $R_{c}$ into account.

We rescale the population density according to

r^{*} = \frac{r}{S^{1 / 3}}, D^{*} (r) = \frac{D (r)}{S^{1 / 3}}

(2)

as proposed in Lemoy and Caruso (2020). It can be seen in Figure 2, that the rescaled curves collapse (Malmgren et al., 2009; Stanley, 1999), i.e. they reasonably well fall onto each other.

Figure 2.

Population density gradients. The rescaled population density is plotted as a function of the rescaled distance to the center, both according to equation (2), for (a) $γ = 2.25$ , (b) $γ = 2.5$ , and (c) $γ = 2.75$ . All panels are in log-log scale. Due to rescaling the values of all realizations and iterations fall onto each other. The shades of blue represent densities, the red line corresponds to the data of an individual curve, and the straight green line is a guide to the eye with slope given by equation (5). The population density asymptotically decays as power-law.

Specifically, we find that the population density decays following a power-law

\frac{D (r)}{S^{1 / 3}} \sim {(\frac{r}{S^{1 / 3}})}^{- α}

(3)

for

r^{*} > r_{p}^{*}

, where

r_{p}^{*}

is the rescaled radius at which the plateau ends and the power-law decay begins. Rescaling does not affect the power-law relation and we conclude that the population density generated by our model follows

D (r) \sim r^{- α}

(4)

for

r^{*} > r_{p}^{*}

. The power-law reasonably well agrees with the empirical results (Lemoy and Caruso, 2020). Since Lemoy and Caruso (2017, 2020) study profiles across many cities at the same time step and we rescale various realizations but across time (instead of only across samples), we hypothesize ergodicity, in a sense that cross-sectional and temporal behavior are the same.

As can be seen in Figure 2, the density gradient exponent α depends on the gravity exponent γ. We find

α = 2 γ - 3

(5)

Small γ-values lead to scattered/sprawled structures and large γ-values lead to compact patches. The value $γ ≃ 2.5$ as estimated for Paris (Rybski et al., 2013) agrees well with $α ≃ 2$ as indicated in Lemoy and Caruso (2020).

Urban fraction gradient

Analogously to the population density, the urban fraction is given by $u (r) = \sum_{k} θ (w_{k}) / \sum_{k} 1$ , where $θ (w_{k}) = 0$ for w_k = 0 and $θ (w_{k}) = 1$ for $w_{k} > 0$ . Again, as proposed in Lemoy and Caruso (2020), we rescale the urban fraction according to

r^{*} = \frac{r}{S^{1 / 2}}, u^{*} (r) = u (r) .

(6)

Similar to the population density, the rescaled curves of urban fraction collapse onto each other in Figure 3. For the urban fraction we find an exponential decay

u (r) = \exp (- \frac{b}{S^{1 / 2}} r)

(7)

for

r^{*} > r_{p}^{*}

. The urban fraction gradient parameter depends on the gravity exponent, i.e.

b / S^{1 / 2} \sim \exp (c γ)

with

c \approx 10 / 3

, see Figure 3(b), (d) and (f).

Figure 3.

Urban fraction gradients. The urban fraction is plotted as a function of the rescaled distance to the center according to equation (6), for (a) + (b) $γ = 2.25$ , (c) + (d) $γ = 2.5$ , and (e) + (f) $γ = 2.75$ . Panels (a), (c), (d) are in lin-lin scale and in panels (b), (d), (f) the vertical axis is logarithmic. Due to rescaling the values of all realizations and iterations fall onto each other. The shades of blue represent densities. The red line in (b), (d), and (f) corresponds to the data of an individual curve, and the straight green line is a guide to the eye from equation (7). The urban fraction decays exponentially.

Equation (7) seems to hold reasonably well (Makse et al., 1995), but overall Lemoy and Caruso (2020) find a slower than exponential decay.

Urban fractality

Next we want to argue that equation (4) is related to fractality (Batty and Longley, 1994; Encarnação et al., 2012; Frankhauser, 2008; Zhou et al., 2017). The fractal dimension d is commonly defined by $M \sim L^{d}$ , i.e. by the way how the mass M of the considered structure changes with linear size L, (see e.g. Bunde and Havlin, 1995). In our case the relation between M and L can be expressed as a mass-radius relation (Daqing et al., 2011; Makse et al., 1998). Moreover, we are studying the density $D = M / L^{2}$ . In combination we can write $D (r) \sim r^{d_{rad} - 2}$ (Batty and Longley, 1994, equation (8.12)). Comparison with equation (4) leads to $α = 2 - d_{rad}$ and with equation (5) to

d_{rad} = 5 - 2 γ

(8)

For $γ = 2.5$ , the resulting structures are zero-dimensional in terms of fractal geometry, i.e. essentially corresponding to a point. For $γ > 2.5$ , we obtain negative fractal dimensions, from which we infer that the mass-radius relation is not valid anymore, i.e. radial symmetry is lost.

Area-perimeter relation

While so far we have studied the resulting w-values of the whole system, from now on we focus on the properties of the major central cluster. To be more specific, here we consider its shape. As introduced by Lovejoy (1982) we first investigate the area-perimeter relation (see also Batty and Longley, 1994: Ch. 6.2), according to which the area $A_{\cdot}$ and the perimeter C of the object under consideration are related by a power-law

C \sim A_{\cdot}^{d_{apr}^{°} / 2}

(9)

where

d_{apr}^{°}

is the fractal dimension of the perimeter. By

A_{\cdot}

we denote the area of the cluster where we fill any empty cells (holes) within the perimeter, accordingly,

A_{c} \leq A_{\cdot}

Figure 4(a) shows an example of the correlations between area and perimeter. As expected there is a power-law relation. We have fitted the exponent in equation (9) based on the evolution of the major central cluster of each realization separately. In Figure 4(b), the resulting fractal dimensions of all realizations are plotted as a function of the various γ-values. There is considerable spreading among the realizations but a minor increase of the average values can be observed from $d_{apr}^{°} \approx 1.25 \dots 1.32$ for $γ = 2.1 \dots 2.75$ , respectively. This range is comparable to the range achieved by the correlated percolation model (CPM) (Makse et al., 1998). Any size dependence of $d_{apr}^{°}$ cannot be studied since in this case size variation defines the dimension.

Figure 4.

Area-perimeter relation. (a) The area and perimeter of the major central cluster are correlated according to a power-law equation (9), here shown for $γ = 2.5$ . In the panel all realizations have been combined. (b) Fractal dimension of the perimeter $d_{apr}^{°}$ according to the area-perimeter relation as a function of the gravity exponent γ. The red dots represent the averages among the realizations.

Box-counting dimensions

Next we employ box-counting to characterize the structure of the major central cluster. The method consists of counting the number of non-overlapping square-shaped boxes necessary to cover an object (see Bunde and Havlin, 1995, and references therein). By varying the size of the box the dimension is quantified via

N_{bc} \sim L^{d}

(10)

where

N_{bc}

is the number of boxes, L their size and d the dimension. We assess the cluster as a whole as well as the envelope of the cluster and denote the dimensions

d_{bc}^{•}

and

d_{bc}^{°}

, respectively.

In Figure 5, we plot the resulting fractal dimensions as a function of the size of the major central cluster $A_{c}$ . It can be seen that the fractal dimensions tend to increase with $A_{c}$ , which is qualitatively consistent with empirical findings and previous results (Rybski et al., 2013; Shen, 2002; Zhou et al., 2017). The correlations are non-linear and more pronounced for $d_{bc}^{•}$ , i.e. the fractal dimension of the entire cluster correlates better with the size.

Figure 5.

Box-counting fractal dimensions of major central cluster and its envelope. The fractal dimensions are plotted vs. cluster size in (a)-(d) for the envelope, $d_{bc}^{°}$ , in (a)+(b), and the entire cluster, $d_{bc}^{•}$ , in (c)+(d). As examples we use (a) + (c) $γ = 2.25$ and (b) + (d) $γ = 2.75$ . The different colors represent the various realizations. The two fractal dimensions are plotted against each other in (e), where it can be seen that the fractal dimension of the envelope is slightly smaller for the larger γ-value. Overall, the influence of the cluster size is stronger than the gravity exponent. Panel (f) illustrates for an example of the population density gradient (semi-log) how the different fractal dimensions are defined by different features of the simulations. The orange curves represents the major central cluster and the green one the entire system, i.e. including small surrounding clusters.

In Zhou et al. (2017), based on 5000 clusters of urban land-cover in Europe, the fractal dimension of the envelope roughly varies between 1.3 and 1.5, while for the cluster itself it varies between 1.3 and 1.7. From our simulations we obtain $d_{bc}^{°}$ roughly between 1.1 and 1.4 (Figure 5(a) and (b)) and $d_{bc}^{•}$ roughly between 1.5 and 1.9 (Figure 5(c) and (d)). However, Zhou et al. (2017) also report an anisometry of the clusters which could affect the fractal dimension. The influence of γ is small and can be seen in Figure 5(e) where we plot the fractal dimensions against each other. The smaller $γ = 2.25$ leads to slightly larger $d_{bc}^{°}$ . Overall, the dependence on the size is more pronounced than the influence of the gravity exponent γ.

It needs to be noted that while $d_{rad}$ in equation (8) describes the fractality of the entire cluster, including the population (i.e. the third dimension), $d_{apr}^{°}$ in equation (9) and $d_{bc}^{°}$ characterize what in the context of the CPM is called percolation front (Makse et al., 1998), i.e. the fuzziness of the envelope (see Figure 5(f)). The measure $d_{bc}^{•}$ is a combination of both, but does not consider the third dimension.

Fundamental allometry

Schläpfer et al. (2015) find a power-law between the average building height and city size. In our context the building height translates into population density so that their relation corresponds to $S_{c} / A_{c} \sim S_{c}^{ϕ}$ . This power-law, in turn, is associated to the fundamental allometry relating the population and area of cities (Batty and Ferguson, 2011; Fluschnik et al., 2016; Nordbeck, 1971; Rybski, 2016; Rybski et al., 2017; Stewart and Warntz, 1958)

S_{c} \sim A_{c}^{δ}

(11)

via

ϕ = 1 - 1 / δ

. Accordingly, in the following we analyze equation (11) for our model, i.e. the major central cluster.

In Figure 6(a) one can see that the resulting populations and areas follow power-laws according to equation (11). The allometry exponent δ depends on the gravity exponent γ, approximately following a parabolic relationship, see Figure 6(b). Schläpfer et al. (2015) report $ϕ ≃ 0.34$ – considering buildings within 2 km from the city center – which corresponds to $δ ≃ 1.52$ and $γ \approx 3$ according to our numerical results. Our model seems to be restricted to $δ > 1$ , which is consistent with the results for the majority of real-world cities (Batty and Ferguson, 2011; Bettencourt and Lobo, 2016).

Figure 6.

Fundamental allometry. The population of the major central cluster as a function of its area is plotted in panel (a) for various γ-values as indicated in the legend. The dotted line has slope 1. Values of all realizations and iterations are shown. The scaling exponent δ according to equation (11) is plotted in panel (b) together with a parabolic regression. The 10 realizations for each γ-value are represented by a box-plot. Larger values of γ lead to increased population density in big cities.

Economics reasoning

We want to motivate equation (4), i.e. propose a setting under which the density gradient goes as

D (r) \sim r^{- α}

(12)

whereas

α \approx 2

. The population density is given by population per area. The area of concentric rings is proportional to the distance from the center, i.e.

\sim 2 π r

. Then

D (r) = \frac{S (r)}{A (r)} \sim \frac{S (r)}{r}

and the population follows

S (r) \sim r^{- α + 1}

. If we normalize to the total population, then we have a probability density

p (r) \sim r^{- α + 1}

(13)

We begin with the common approach according to which the commodity Z is given by the wages W minus the housing rent R and the transportation costs T (see e.g. Barthelemy, 2016: Ch. 3.3), i.e.

Z = W - R - T

(14)

which is maximized by minimizing expenses

max (Z) = W - min (R + T)

(15)

In the mono-centric case, it is common sense that the rents decrease further away from the city center but the transportation costs increase. We assume power-law relations

R = a r^{- ρ}

(16)

T = b r^{τ}

(17)

with the distance r from the center and

ρ, τ > 0

. The power-law rent profile is motivated below, for the power-law transportation costs (see e.g. Fabinger, 2012: Sec. 1.4.3) and references therein. Delloye et al. (2018) use linear transportation costs (τ = 1).

In order to find the optimal distance to the center, the derivative of R + T needs to be zero, i.e.

r_{opt} = {(\frac{ρ a}{τ b})}^{1 / (ρ + τ)}

(18)

By summing R with T we are essentially comparing apples with oranges. However, the prefactors a and b determine the weights they have relative to each other. Certainly, for wealthy people, the rent becomes less of an issue and the weight should be smaller while transportation is similar for everyone (people spend 20% to 30% of their time commuting (Kahneman et al., 2006)). We assume $a \sim W^{- 1}, b = const$ (see Figure 7) and obtain

r_{opt} \sim W^{- 1 / (ρ + τ)}

(19)

Figure 7.

Illustration of the sum of housing rent and transportation costs as well as the influence of the weight a. The values ρ = 2, τ = 3, and b = 1 have been used exemplary. With decreasing a the minimum moves towards the center.

Wealthier people can afford living close to the city center while low income population is pushed outward.

Further, we take the power-law income or wealth distribution

p (W) \sim W^{- ζ}

(20)

with

ζ \approx 2.5

(

ζ_{USA} ≃ 2.4

(Brzezinski, 2014; Levy and Solomon, 1997)). The exponent ζ is also related to the Gini coefficient via

ζ = \frac{1}{2 G} + \frac{3}{2}

(Pfähler, 1985). Typical values are between G = 0.65 and G = 0.80 corresponding to

ζ ≃ 2.27

and

ζ ≃ 2.13

, respectively.

If two quantities A and B follow power-law distributions with pdfs $p (A) \sim A^{- ζ_{A}}$ and $p (B) \sim B^{- ζ_{B}}$ , then the transformation $B \sim A^{β}$ with $β = (ζ_{A} - 1) / (ζ_{B} - 1)$ translates one into the other (Gomez-Lievano et al., 2012). Comparison leads to $ζ_{A} = ζ, ζ_{B} = α - 1$ , and

\frac{- 1}{ρ + τ} = \frac{ζ - 1}{α - 2}

(21)

The lhs is negative and since $ζ ≫ 1$ the rhs can only become negative if $α < 2$ . At α = 2 a transition occurs and for $α > 2$ the mono-centric assumption does not hold. This is consistent with the transition at $γ = 2.5$ in equation (8). From $α = 2 - (ζ - 1) (ρ + τ)$ it can be seen that only small values of τ and ρ lead to α close to 2.

As a critical remark, we need to add that it is not a surprise to obtain a power-law (or a relation between exponents) when the derivation itself is based on power-laws. However, in economics power-laws are theoretically understood and empirically established (Gabaix, 2016).

Housing rent

Here we want to motivate equation (16) and the exponent ρ. The housing rent is determined by supply and demand. We postulate that the number of people willing to pay rent larger than R decreases as a power-law with R

P_{D} (X \geq R) ≃ R^{- ϵ_{D}}

(22)

Analogously, the number of people willing to sell property or rent it out for a price lower than R decreases as a power-law with R

P_{S} (X \leq R) ≃ 1 - R^{- ϵ_{S}}

(23)

The market price is then given by the price where both curves cross each other

a_{D} P_{D} (X \geq R_{\times}) = a_{S} P_{S} (X \leq R_{\times})

(24)

a_{D} R_{\times}^{- ϵ_{D}} + a_{S} R_{\times}^{- ϵ_{S}} = a_{S}

(25)

where the factors

a_{D}

and

a_{S}

are required to adjust for the amount and convert the cdfs into cumulative frequency distributions. Increasing availability should decrease the price and increasing demand should increase the price. We assume

a_{D} \sim 1 / r

and

a_{S} \sim r

leading to

1 / r R_{\times}^{- ϵ_{D}} + r R_{\times}^{- ϵ_{S}} \sim r

. For large r the second term dominates, implying

R_{\times} = const

. Thus, large supply leads to a (low) price that is independent of the location. For small r the first term dominates, leading to

R_{\times} \sim r^{- 2 / ϵ_{D}}

(26)

i.e. $ρ = 2 / ϵ_{D}$ . The price is dominated by the demand.

Linearity should work if we consider the area. In case of living space/apartments another exponent might be necessary to take changes of density into account, i.e. $a_{D} \sim r^{- δ_{D}}$ and $a_{S} \sim r^{δ_{S}}$ . With the same reasoning as before, we then obtain $R_{\times} \sim r^{- (δ_{S} + δ_{D}) / ϵ_{D}}$ . In particular, if we consider equation (13) and $δ_{D} = δ_{S} = α - 1$ , then $ρ = (2 α - 2) / ϵ_{D}$ . For $α \approx 2$ we have $ρ \approx 2 / ϵ_{D}$ .

Summary and discussion

In summary, our simulations show that the gravitational approach – according to which the probability of incremental growth is proportional to $d^{- γ}$ – is capable of reproducing radial gradients of real-world cities. We numerically find a relation between the gravity exponent γ and the population density exponent α, suggesting equivalence i.e. the power-law population gradient is an expression of the gravitation (or vice versa). Accordingly, our results confirm the idea of a friction of distance (Batty and Ferguson, 2011; Benenson and Torrens, 2004: Sec. 3.2.2). However, the strength of proximity follows a power-law and it can be anticipated that an exponential functions instead of $d^{- γ}$ in equation (1) will not lead to rich spatial complexity. If we accept that the model generates structures that resemble real-world cities, then we can conclude that gravitation represents a composite mechanism of the various attractive processes influencing location choices (proximity to friends and work, availability of infrastructure, clustering of business types, etc.).

Our results add to the pioneering work by Batty and Sik Kim (1992) who described a power-law population density gradient – in contrast to an exponential one (Clark, 1951). However, the proposed range of α between 0 and 1 corresponds to γ between 1.5 and 2, which is below the range investigated here. For $γ < 2$ the emerging structures are too noisy (Rybski et al., 2013) and unrealistic compared to real-world cities. A recently suggested value is $α ≃ 2$ (Lemoy and Caruso, 2017) which corresponds to $γ ≃ 2.5$ (Rybski et al., 2013). Interestingly, at this precise value of α = 2 and $γ = 2.5$ the fractal dimension is $d_{rad} = 0$ which agrees with the perception of cities as (zero-dimensional) points. For $γ > 2.5$ a transition occurs where the fractal dimension is not defined or does not well-behave which we denote singularity.

Moreover, it needs to be mentioned that we study our model results in terms of mono-centric cities. If the main cluster merges with surrounding smaller ones, then sub-centers can appear, but overall the main center dominates, as illustrated in Figure 1. It remains to be studied what happens in the regime $γ > 2.5$ . In any case, a coherent definition and an appropriate measure for poly-centrism are lacking.

We investigate additional three fractal dimensions characterizing the structure of the major central cluster, disregarding the population density. Overall we find that the fractality is dominated by the size of the cluster while the gravity exponent γ has a minor influence. This is consistent with various previous papers.

Our approach also leads to urban allometry between population and area, although the scaling seems to be restricted to $δ > 1$ in equation (11), i.e. the case where large cities exhibit higher densities. In a sense the exponent γ determines how sprawled or compact the emergent cities are. If one could find a policy instrument to control γ, then one could use it to influence the urban development in the desired way. This could address sustainability questions, e.g. related to the ratio of land consumption rate to population growth. Specifically, larger γ-values lead to more compact cities – due to the fundamental allometry equation (11) this influences mostly the large ones.

An alternative model that elegantly generates spatial complexity and radial gradients is diffusion-limited aggregation (DLA) (Batty, 2013; Batty and Longley, 1994; Batty et al., 1989; Fotheringham et al., 1989; Witten and Sander, 1981). Contrasting equation (7), DLA leads to a power-law gradient of the urban fraction (Fotheringham et al., 1989, equation (3)). A form of allometry, equation (11), is also obtained from DLA (Fotheringham et al., 1989, equation (5)). The fractal dimension of the DLA in its basic form is $\approx 1.71$ (e.g. Batty and Longley, 1994). However, as the present model also grows in the third dimension, DLA can rather be compared to the binary gravitational model (Rybski et al., 2013).

In contrast to the correlated percolation model (CPM) proposed in (Makse et al., 1995, 1998), where the urban fraction gradient and the structure are introduced artificially, in the gravitational approach presented here they are emergent. Moreover, it is not straight forward to extend the CPM to also simulate population density. It would be interesting to analyze which gradients are generated by the spatial network model (SNM) (Frasco et al., 2014; Wickramasinghe et al., 2018).

The q_j in equation (1) are often interpreted as potential of a gravitational force $F_{j, k}$ (Batten and Boyce, 1987). In physics, they are related via $F = - \nabla q$ and as the distance appears $\sim d^{- γ}$ in the potential, it should be $\sim d^{- (γ + 1)}$ in the force. However, the analogy only works partly. First, the system does not have any dynamics and the kinetic energy as counterpart to the potential is missing. Second, the q_j are probabilities and more similar to the probability density of finding a particle at a given place, i.e. the squared modulus of a wave function, but it is questionable if a wave function makes sense in this context.

We also would like to discuss some limitations of our gravitational model. (i) The maximum urban fraction reaches 1 which is higher than in real cities. Analogously, population shows unbounded growth in the core and not a plateau as in real cities (although in Figure 2 a plateau can be seen in the log-log scale, in lin-lin representation, it is negligible). (ii) The assumed proportionalities between the w-values, population density and building height do not affect the model interpretations but for the comparison with real-world cities they represent rough assumptions and might require refinements from follow-up studies (Biljecki et al., 2016). In particular, a central business district and similar features would require to distinguish residential from commercial and other uses. (iii) Real-world cities are rarely radial and many exhibit anisometry (Zhou et al., 2017), which in most cases results from landscape constraints (e.g. coast lines). Our model apparently does not reproduce such anisometry but was also not intended to do so. (iv) The growth is exogenous and the constant growth parameter g leads to idealized urban development trajectories.

In principle the model can be extended by another exponent ϵ, i.e. $w_{k}^{ϵ} d_{j, k}^{- γ}$ in equation (1), giving more or less dense cells more or less weight. For the sake of simplicity we did not follow this approach. Moreover, in the context of complex networks it has been shown that ‘nonlinear preferential attachment’ ( $ϵ \neq 1$ ) leads to degree distributions which are different from power-laws (Krapivsky et al., 2000). For systems of cities this would imply deviations from Zipf ’s law.

Last but not least, we would like to discuss an outlook to future work. (i) Recently, Volpati and Barthelemy (2018) proposed a dispersion index to characterize the degree of localization in populous areas. A direction of future research could be to apply it to our model output and establish a relation. (ii) It could also be of interest to operationalize the gravitational approach in order to apply it to real-world data (Jones and O’Neill, 2016). More landscape features need to be taken into account for a realistic modeling. (iii) The described gradients could be related to other quantities, such as the Urban Heat Island (UHI) effect (Watkins et al., 2002: Figure 6; Zhou et al., 2017).

Footnotes

Acknowledgements

We would like to thank A. Brenner, Y. Liu, S. Thies, and B. Zhou for useful discussions. This work emerged from ideas discussed at the symposium Cities as Complex Systems (Hanover, July 13th-15th, 2016) which was generously founded by Volkswagen Foundation.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Y. Li thanks China Scholarship Council (CSC) for financial support. D. Rybski is grateful to the Leibniz Association (project IMPETUS) for financially supporting this research.

Yunfei Li is a PhD student of Urban Transformations group at Potsdam Institute for Climate Impact Research (PIK). He majored in geographic information system during undergraduate and in water conservancy information technology for Masters. Now he is working in the field of urban climate modelling to study the impact of urbanization on urban climate and the interaction between urban heat island effect and urbanization under the background of climate change.

Diego Rybski is a senior scientist and the deputy head of the Urban Transformations group at the Potsdam Institute for Climate Impact Research (PIK). At University of Potsdam he lectures on complex systems and sustainability. He wrote his dissertation about applying methods from statistical physics to climatological time series. After two years of postdoctoral research at the City College of New York, in recent years he dedicates his research to the understanding of cities and their role in the climate change discourse (mitigation&adaptation).

Juergen P. Kropp studied chemistry and physics and holds a PhD in theoretical physics. He heads the Urban Transformations group at Potsdam Institute for Climate Impact Research and acts as deputy chair for Research Domain II “Climate Resilience”. At the University of Potsdam he holds an extraordinary professorship for Climate Change and Sustainable Development. He coordinated many international projects on climate change adaptation and mitigation in cities. Juergen's research interests are devoted to complex systems analysis in relation to climate change sustainability challenges.

References

Barthelemy

(2016) The Structure and Dynamics of Cities – Urban Data Analysis and Theoretical Modeling. Cambridge: Cambridge University Press.

Batten

Boyce

(1987) Spatial interaction, transportation, and interregional commodity flow models. Handbook of Regional and Urban Economics 1: 357–406.

Batty

(2013) The New Science of Cities. Cambridge: MIT Press.

Batty

Ferguson

(2011) Defining city size. Environment and Planning B 38(5): 753–756.

Batty

Longley

(1994) Fractal Cities: A Geometry of Form and Function. San Diego: Academic Press Inc.

Batty

Longley

Fotheringham

(1989) Urban growth and form: Scaling, fractal geometry, and diffusion-limited aggregation. Environment and Planning A 21(11): 1447–1472.

Batty

Sik Kim

(1992) Form follows function: Reformulating urban population density functions. Urban Studies 29(7): 1043–1069.

Benenson

Torrens

(2004) Geosimulation: Automata-Based Modeling of Urban Phenomena. West Sussex: John Wiley & Sons Ltd.

Bettencourt

LMA

Lobo

(2016) Urban scaling in Europe. Journal of the Royal Society Interface 13(116): 20160005.

10.

Biljecki

Ohori

Ledoux

; et al. (2016) Population estimation using a 3d city model: A multi-scale country-wide study in the Netherlands. PLoS One 11(6): e0156808.

11.

Brzezinski

(2014) Do wealth distributions follow power laws? Evidence from ‘rich lists’. Journal of Physics A 406: 155–162.

12.

Bunde

Havlin

(1995) Chapter 1. In: Fractals in Science. A Brief Introduction to Fractal Geometry. Berlin: Springer-Verlag, pp. 1–25.

13.

Carrothers

GAP

(1956) An historical review of the gravity and potential concepts of human interaction. Journal of the American Institute of Planners 22(2): 94–102.

14.

Clark

(1951) Urban population densities. Journal of the Royal Statistical Society 114(4): 490–496.

15.

Daqing

Kosmidis

Bunde

; et al. (2011) Dimension of spatially embedded networks. Nature Physics 7(6): 481.

16.

Delloye

Lemoy

Geoffrey

(2018) Alonso and the scaling of urban profiles. arXiv.org e-Print archive arXiv:1801.07512 [physics.soc-ph]. Available at: https://arxiv.org/abs/1801.07512 (accessed 5 April 2019).

17.

Encarnação

Gaudiano

Santos

; et al. (2012) Fractal cartography of urban areas. Scientific Reports 2(527): srep00527.

18.

Fabinger

(2012) Essays on trade and imperfectly competitive markets. Doctoral Dissertation, Harvard University. Available at: http://nrs.harvard.edu/urn-3:HUL.InstRepos:9548614 (accessed 4 April 2019).

19.

Fluschnik

Kriewald

Ros

AGC

; et al. (2016) The size distribution, scaling properties and spatial organization of urban clusters: A global and regional percolation perspective. International Journal of Geo-Information 5(7): 110.

20.

Fotheringham

Batty

Longley

(1989) Diffusion-limited aggregation and the fractal nature of urban growth. Papers in Regional Science Association 67(1): 55–69.

21.

Frankhauser

(2008) Fractal geometry for measuring and modelling urban patterns. In: Albeverio

Andrey

Giordano

; et al. (eds) The Dynamics of Complex Urban Systems. Heidelberg: Springer, pp. 213–243.

22.

Frasco

Sun

Rozenfeld

; et al. (2014) Spatially distributed social complex networks. Physical Review X 4(1): 011008.

23.

Gabaix

(2016) Power laws in economics: An introduction. Journal of Economic Perspectives 30(1): 185–206.

24.

Gomez-Lievano

Youn

Bettencourt

LMA

(2012) The statistics of urban scaling and their connection to Zipf’s law. PLoS One 7(7): e40393.

25.

Guérois

Pumain

(2008) Built-up encroachment and the urban field: A comparison of forty European cities. Environment and Planning B 40(9): 2186–2203.

26.

Hoshen

Kopelman

(1976) Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm. Physical Review B 14(8): 3438–3445.

27.

Jones

O’Neill

(2016) Spatially explicit global population scenarios consistent with the shared socioeconomic pathways. Environmental Research Letters 11(8): 084003.

28.

Kahneman

Krueger

Schkade

; et al. (2006) Would you be happier if you were richer? A focusing illusion. Science 312(5782): 1908–1910.

29.

Krapivsky

Redner

Leyvraz

(2000) Connectivity of growing random networks. Physical Review Letters 85(21): 4629.

30.

Kriewald

Fluschnik

Reusser

; et al. (2016) This is a test entry of type @ONLINE. Available at: https://cran.r-project.org/web/packages/osc/index.html (accessed 5 April 2019).

31.

Lemoy

Caruso

(2017) Scaling evidence of the homothetic nature of cities. arXiv.org e-Print archive arXiv:1704.06508 [physics.soc-ph]. Available at: https://arxiv.org/abs/1704.06508 (accessed 5 April 2019).

32.

Lemoy

and Caruso

(2020) Evidence for the homothetic scaling of urban forms. Environment and Planning B: Urban Analytics and City Science 47(5): 870–888.

33.

Levy

Solomon

(1997) New evidence for the power-law distribution of wealth. Journal of Physics A 242(1-2): 90–94.

34.

Lovejoy

(1982) Area-perimeter relation for rain and cloud areas. Science 216(4542): 185–187.

35.

Makse

Andrade

Batty

; et al. (1998) Modeling urban growth patterns with correlated percolation. Physical Review E 58(6): 7054–7062.

36.

Makse

Havlin

Stanley

(1995) Modeling urban-growth patterns. Nature 377(6550): 608–612.

37.

Malmgren

Stouffer

Campanharo

ASLO

; et al. (2009) On universality in human correspondence activity. Science 325(5948): 1696–1700.

38.

Morgan

(2004) The exaggerated death of geography: Learning, proximity and territorial innovation systems. Journal of Economic Geography 4(1): 3–21.

39.

Nordbeck

(1971) Urban allometric growth. Geografiska Annaler: Series B 53(1): 54–67.

40.

Peiravian

Kermanshah

Derrible

(2014) Spatial data analysis of complex urban systems. In: 2014 IEEE international conference on big data. Piscataway: IEEE, pp. 1–6.

41.

Pfähler

(1985) Relative concentration curve: Functional form and measures of non-proportionality. Bulletin of Economic Research 37(3): 201–211.

42.

Rozenfeld

Rybski

Andrade

Jr ; et al. (2008) Laws of population growth. In: Proceedings of the national academy of sciences of the United States of America 105(48): 18702–18707.

43.

Rozenfeld

Rybski

Gabaix

; et al. (2011) The area and population of cities: New insights from a different perspective on cities. American Economic Review 101(5): 2205–2225.

44.

Rybski

(2016) Relating urban scaling, fundamental allometry, and density scaling. arXiv.org e-Print archive arXiv:1609.01217 [physics.soc-ph]. Available at: https://arxiv.org/abs/1609.01217 (accessed 5 April 2019).

45.

Rybski

Reusser

Winz

; et al. (2017) Cities as nuclei of sustainability? Environment and Planning B 44(3): 425–440.

46.

Rybski

Ros

AGC

Kropp

(2013) Distance weighted city growth. Physical Review E 87(4): 042114.

47.

Schläpfer

Lee

Bettencourt

LMA

(2015) Urban skylines: Building heights and shapes as measures of city size. arXiv.org e-Print archive arXiv:1512.00946 [physics.soc-ph]. Available at: https://arxiv.org/abs/1512.00946 (accessed 5 April 2019).

48.

Shen

(2002) Fractal dimension and fractal growth of urbanized areas. International Journal of Geographical Information Science 16(5): 419–437.

49.

Stanley

(1999) Scaling, universality, and renormalization: Three pillars of modern critical phenomena. Reviews of Modern Physics 71(2): S358–S366.

50.

Stewart

(1948) Demographic gravitation: Evidence and applications. Sociometry 11(1/2): 31–58.

51.

Stewart

Warntz

(1958) Physics of population distribution. Journal of Regional Science 1(1): 99–121.

52.

Volpati

Barthelemy

(2018) The spatial organization of the population density in cities. arXiv.org e-Print archive arXiv:1804.00855 [physics.soc-ph]. Available at: https://arxiv.org/abs/1804.00855 (accessed 5 April 2019).

53.

Watkins

Palmer

Kolokotroni

; et al. (2002) The London heat island: Results from summertime monitoring. Building Services Engineering Research and Technology 23(2): 97–106.

54.

White

Engelen

Uljee

(2015) Modeling Cities and Regions as Complex Systems: From Theory to Planning Applications. Cambridge: MIT Press.

55.

Wickramasinghe

Onyerikwu

Sun

; et al. (2018) Modeling spatial social complex networks for dynamical processes. Complexity 2018: 1428719.

56.

Witten

Jr Sander

(1981) Diffusion-limited aggregation, a kinetic critical phenomenon. Physical Review Letters 47(19): 1400–1403.

57.

Zhou

Rybski

Kropp

(2017) The role of city size and urban form in the surface urban heat island. Scientific Reports 7(1): 4791.

58.

Zipf

(1946) The

p_{1} p_{2} / d

hypothesis: On the intercity movement of persons. American Sociological Review 11(6): 677–686.