Effects of changing population or density on urban carbon dioxide emissions

TL;DR

This paper introduces a comprehensive framework to analyze how population and area, along with their interactions, influence urban CO2 emissions, revealing that population changes have a greater impact than density variations.

Contribution

It develops a generalized model considering population, area, and their interactions, improving emission predictions and uncovering the coupled effects of urban metrics.

Findings

Population changes have a greater effect on emissions than density changes.

Larger cities show a higher impact of population or density variations on emissions.

The model accounts for initial values of population and area in emission variations.

Abstract

The question of whether urbanization contributes to increasing carbon dioxide emissions has been mainly investigated via scaling relationships with population or population density. However, these approaches overlook the correlations between population and area, and ignore possible interactions between these quantities. Here, we propose a generalized framework that simultaneously considers the effects of population and area along with possible interactions between these urban metrics. Our results significantly improve the description of emissions and reveal the coupled role between population and density on emissions. These models show that variations in emissions associated with proportionate changes in population or density may not only depend on the magnitude of these changes but also on the initial values of these quantities. For US areas, the larger the city, the higher is the impact of changing its population or density on its emissions; but population changes always have a greater effect on emissions than population density.

Tables

Table 1. Supplementary Table 1: Exponents associated with the urban scaling ( β 𝛽 \beta ) and per capita density scaling ( α 𝛼 \alpha ). Abbreviations: Reduced Major Axis Regression (RMA), Ordinary Least Square Regression (OLS), Prais-Winsten Regression (PW).

Method, reference	Country	Year	Urban Units	Sample Size	Sector	Value (Lower,Upper)
						$\beta$
RMA 20	Annex I	2005	Large Cities	22	N/A	$0.87(0.59,0.95)$
RMA 20	Non- Annex I	2005	Large Cities	39	N/A	$1.18(0.90,1.49)$
OLS 19	US	1999- 2011	Metropolitan Areas	1875	Local CO2 and five air pollutants	$0.75(0.71,0.79)$
OLS 16	US	1982- 2011	Urban Centers	101	Excess CO2 from congestion	$1.27(1.21,1.33)$
PW 16	US	1982- 2011	Urban Centers	101	Excess CO2 from congestion	$1.14(1.38,1.44)$
OLS 17	US	N/A	N/A	122	N/A	$0.92(0.90,0.95)$
OLS 17	India	N/A	N/A	42	N/A	$1.13(1.02,1.23)$
OLS 17	UK	N/A	N/A	35	N/A	$0.76(0.68,0.83)$
OLS 13	GB	2010	City & Administrative Boundaries	28	Transport CO2 emissions	$1.02(0.96,1.08)$
OLS 11	US	2010	Metropolitan Areas	101	Excess CO2 from congestion	$1.26(1.17,1.35)$
OLS 11	OECD Countries	2014	Metropolitan Areas	268	Transport CO2 emissions	$1.21(1.11,1.31)$
OLS 9	US	2002	CCA clusters ($l=5$ km, $D=1000$)	2281	All Vulcan sectors	$1.38(1.35,1.41)$
OLS 9	US	2002	CCA clusters ($l>10$ km, )	N/A	All Vulcan sectors	$1.46(1.44,1.48)$
OLS 8	US	1999- 2008	Metropolitan Areas	933	All Vulcan sectors	$0.93(0.92,0.95)$
						$\alpha$
OLS 15	US	2000	CCA/GRUMP/ GLC ($l=1$ km)	4585	Buildings and on-road emissions	$-0.78(\text{N/A},\text{N/A})$
OLS 15	US	2000	CCA/GRUMP/ GLC ($l=5$ km)	3285	Buildings and on-road emissions	$-0.79(\text{N/A},\text{N/A})$
OLS 15	US	2000	CCA/GRUMP/ GLC ($l=10$ km)	2786	Buildings and on-road emissions	$-0.82(\text{N/A},\text{N/A})$
OLS 15	US	2000	CCA/GRUMP ($l=1$ km)	5182	Buildings and on-road emissions	$-0.90(\text{N/A},\text{N/A})$
OLS 15	US	2000	CCA/GRUMP ($l=5$ km)	2156	Buildings and on-road emissions	$-1.10(\text{N/A},\text{N/A})$
OLS 15	US	2000	CCA/GRUMP ($l=10$ km)	1538	Buildings and on-road emissions	$-1.13(\text{N/A},\text{N/A})$
OLS 5	Global	1980	World cities	32	Gasoline consumption	$-0.92(-0.85,-0.99)$
OLS 23	EU	2007-2009	European cities	134	Transport CO2 emissions	$-0.19(\text{N/A},\text{N/A})$
OLS 23	EU	2007-2009	France	24	Transport CO2 emissions	$-0.39(\text{N/A},\text{N/A})$
OLS 23	EU	2007-2009	Spain	22	Transport CO2 emissions	$-0.40(\text{N/A},\text{N/A})$

Table 2. Supplementary Table 2: Summary of the main properties of Cobb-Douglas [Eq. ( 3 )] and translog [Eq. ( 5 )] models. All marginal products are expressed in logarithmic scale.

Property	Cobb-Douglas model, Eq. (3)	Translog model, Eq. (5)
Elasticity of scale in terms of $P$ and $A$	${\beta }_P+{\beta }_A$	${\beta }_P+{\beta }_A+{\beta }_C\log \left(PA\right)$
Elasticity of scale in terms of $P$ and $P/A$	${\beta }_P$	${\beta }_P+{\beta }_C\log A$
Technical rate of substitution between $P$ and $A$	$-\frac{{\beta }_P}{{\beta }_A\left(P/A\right)}$	$\frac{-1}{\left(P/A\right)}\left(\frac{{\beta }_P+{\beta }_C\log A}{{\beta }_A+{\beta }_C\log P}\right)$
Elasticity of substitution	1	$-\frac{{\beta }_P+{\beta }_A+{\beta }_C\log \left(PA\right)}{C\left({\beta }_P+{\beta }_A+{\beta }_C\log \left(PA\right)-2{\beta }_C\right)}$
Marginal product of population (in terms of $P$ and $A$)	${\beta }_P$	${\beta }_P+{\beta }_C\log A$
Marginal product of area (in terms of $P$ and $A$)	${\beta }_A$	${\beta }_A+{\beta }_C\log P$
Marginal product of population (in terms of $P$ and $P/A$)	${\beta }_P+{\beta }_A$	${\beta }_P+{\beta }_A+{\beta }_C\log \left(PA\right)$
Marginal product of density (in terms of $P$ and $P/A$)	$-{\beta }_A$	$-{\beta }_A-{\beta }_C\log P$
Condition for population dominates over density	$\left|{\beta }_P+{\beta }_A\right|>\left|{\beta }_A\right|$	$A>{10}^{-{\beta }_P/{\beta }_A}$ (if ${\beta }_P,{\beta }_A,{\beta }_C>0$ and $P,A>1$)
Condition for decreasing returns to scale (in terms of $P$ and $A$)
Condition for decreasing returns to scale (in terms of $P$ and $P/A$)

Equations

$$C\sim P^{\beta },$$

$$C/P\sim (P/A{)}^{\alpha },$$

$$C\sim P^{{\beta }_P}{A}^{{\beta }_A}\text{or}\log C\sim {\beta }_P\log P+{\beta }_A\log A,$$

$$C\sim ({\beta }_P{P}^{-\gamma }+{\beta }_A{A}^{-\gamma }{)}^{-1/\gamma },$$

$$\log C\sim {\beta }_P\log P+{\beta }_A\log A+{\beta }_C\log (P)\log (A),$$

$${\Omega }^{\ast }=10^{\frac{1-{\beta }_P-{\beta }_A}{{\beta }_C}}.$$

$$\begin{array}{rl}C& \sim P^{{\beta }_P+{\beta }_A+{\beta }_C\log P}(P/A{)}^{-{\beta }_A-{\beta }_C\log P},\end{array}$$

$$C\sim P^{\beta }=P^{{\beta }_1}{P}^{{\beta }_2},$$

$$C\sim P^{{\beta }_1}{A}^{{\beta }_2\delta },$$

$$\beta ={\beta }_P+\frac{{\beta }_A}{\delta },$$

$${\beta }_P=\frac{\beta \delta -1}{\delta -1}\text{and}{\beta }_A=\frac{\delta (\beta -1)}{1-\delta },$$

$$C\sim P^{\theta }{A}^{1-\theta },$$

$$C/A\sim (P/A{)}^{\theta },$$

$$y\sim a_1{x}_1+a_2{x}_2,$$

$$\log C\sim {\tilde{\beta }}_P\widehat{\log P}+{\tilde{\beta }}_A\widehat{\log A},$$

$$\widehat{\log P}=\frac{\log P-{\mu }_P}{{\sigma }_P}\text{and}\widehat{\log A}=\frac{\log A-{\mu }_A}{{\sigma }_A}.$$

$${\beta }_P=\frac{{\tilde{\beta }}_P}{{\sigma }_P{\sigma }_A}\text{and}{\beta }_A=\frac{{\tilde{\beta }}_A}{{\sigma }_P{\sigma }_A},$$

$$\log C\sim {\tilde{\beta }}_P\widehat{\log P}+{\tilde{\beta }}_A\widehat{\log A}+{\tilde{\beta }}_C\widehat{\log (P)}\widehat{\log (A)},$$

$${\beta }_P=\frac{{\sigma }_A{\tilde{\beta }}_P-{\mu }_A{\tilde{\beta }}_C}{{\sigma }_P{\sigma }_A},{\beta }_A=\frac{{\sigma }_P{\tilde{\beta }}_A-{\mu }_P{\tilde{\beta }}_C}{{\sigma }_P{\sigma }_A},\text{and}{\beta }_C=\frac{{\tilde{\beta }}_C}{{\sigma }_P{\sigma }_A},$$

$$\begin{array}{rl}\log C& \sim {\beta }_P\log P+{\beta }_A\log A+{\beta }_C\log (P)\log (A)\\ & +{\beta }_P^{\prime }(\log P{)}^2+{\beta }_A^{\prime }(\log A{)}^2,\end{array}$$

Full text

Introduction

Carbon dioxide (CO2) emissions are considered one of the main causes of Earth’s climate change 1. Despite covering only 0.4 %-0.9 % of global land surfaces 2, urban areas are responsible for more than 70% of such emissions 3; 4. This fact assigns cities a central role in pursuing solutions and mitigation strategies for the global climate change problem. Because of that, researchers from several disciplines have investigated the effects of urbanization on CO2 emissions 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20. The question of whether urbanization promotes or mitigates climate change is ubiquitous among these works, and the approaches to probe such issues differ, but can be roughly organized into two groups.

Underlying the first approach, there is the so-called urban scaling hypothesis 21; 22, which states that city emissions ($C$) are described by a power-law function of population size ($P$), that is, $C\sim P^{\beta}$, where $\beta$ is the urban scaling exponent or the scale-invariant elasticity. For CO2 emissions in US, researchers have reported a 1:1 relationship ($\beta\approx 1$, constant returns to scale) with the population of metropolitan areas 8, while the same quantity was found to scale superlinearly ($\beta=1.46$, increasing returns to scale) when defining the US cities as connected urban spaces 9. When considering local air pollution, an exponent $\beta\approx 3/4$ (decreasing returns to scale) was observed for US metropolitan areas 19. There is also evidence supporting the idea that the scaling between emissions and population depends on the degree of economic development of the urban systems, with increasing returns to scale ($\beta>1$) observed for cities of developing countries and economy of scale ($\beta<1$) for developed ones 17. The second approach is focused on the understanding how population density affects CO2 emissions per capita 5; 6; 23; 12; 14; 15, that is, to investigate the relationship between $C/P$ and $P/A$, where $A$ stands for the urban unit area. A recent work has proposed that CO2 emissions per capita (related to buildings and on-road sectors) and population density are related via the power-law $C/P\sim(P/A)^{\alpha}$, with an exponent $\alpha\approx-0.8$ for the US urban areas 15.

Although these two bodies of the urban CO2 literature are strongly linked by the purpose of understanding how urbanization affects climate change, they have operated widely independent from each other, and their approaches are perceived as different issues. Researchers using urban scaling are assuming population size as the most relevant urban feature for describing CO2 emissions, while those working with the per capita density scaling consider population density as the most significant covariate. Both approaches, however, have produced controversial results regarding the influence of population or population density on urban emissions (see, for instance, Refs. 7; 12 and Supplementary Table 1). Large part of these discrepancies can be attributed to different methodologies for estimating CO2 emissions and defining the boundaries of urban areas, but also because both approaches ignore that population and area are correlated 24; 25 and the influence of a possible interconnected role between these quantities on urban emissions.

Inspired by the economic theory of production functions 26, we propose here a new approach for investigating emissions in urban areas that simultaneously considers the effects of population and area along with possible interactions between these urban metrics. We show that our models recover the two conventional approaches when ignoring the effects of urban area (urban scaling) or when assuming that the emissions display constant returns to scale with population and area (per capita density scaling). When compared with the two conventional approaches, our models provide a significantly better description for the emissions in US urban areas. These results confirm the predictive power of the interactions between population and area, which in turn have intriguing consequences about the effect of these quantities on urban emissions. Our approach indicates that emissions may display decreasing or increasing returns to scale with population and area depending on whether the product $P\times A$ exceeds a particular threshold. We further find that the impact of a proportionate change in the population and density of a city on its emissions increase with its area but always have decreasing returns to scale; moreover, changes in population always have more impact on emissions than changes in density.

Results

Urban scaling and per capita density scaling of emissions

We start by revisiting how population scaling and per capita density scaling approaches have been applied for investigating CO2 emissions in urban areas. To do so, we used the same dataset reported by Gudipudi et al. 15 which comprises CO2 emissions (sum of on-road and building emissions) in US urban areas in the year 2000. As described in Methods, this dataset is constructed by combining gridded data from different sources, and by applying the city clustering algorithm 27 for defining the urban units. There are a total of 3285 urban units and for each one we have population ($P$ in raw counts), area ($A$ in km2), and CO2 emissions ($C$ in tonnes of CO2).

Having defined our variables, within the urban scaling framework, CO2 emissions and population size are related via the power-law function

[TABLE]

where $\beta$ is the scaling exponent. To estimate the parameter $\beta$, we have applied the usual least-squares method to the relationship between $\log C$ and $\log P$. This approach leads to $\beta=0.48\pm 0.01$ ($p$-value $=0$, permutation test, Supplementary Figure 1) and the relationship between both variables (on logarithm scale) is shown in Fig. 1A. At a cursory glance, this value of $\beta$ indicates a sublinear trend between emissions and population, so that a 1% increase in the population level of a city associates with only 0.48% increase in its CO2 emissions. However, a closer inspection of Fig. 1A shows that Eq. (1) deviates systemically from the data and underestimates the emissions in large populated areas. In addition to that, the exponent $\beta$ is likely to be biased by the confounding effect of area because population and area of urban units are correlated to each other via a power-law relation 24; 25.

On the other hand, within the per capita density scaling framework, the relationship between CO2 emissions per capita and population density is described by the power-law function 15

[TABLE]

where $\alpha$ is another scaling exponent. Figure 1B illustrates this relationship on logarithmic scale ($\log C/P$ versus $\log P/A$) from which we have estimated $\alpha=-0.79\pm 0.01$ ($p$-value $=0$, permutation test, Supplementary Figure 1) via ordinary-least-squares method. We observe that Eq. (2) is slightly better than Eq. (1) for describing our data (Supplementary Figure 2), and does not seem to exhibit any systematic deviation. Our estimate for $\alpha$ is in agreement with the results reported by Gudipudi et al. 15, and indicates that every $1\%$ increase in the population density of a city associates with a $0.79\%$ reduction in its CO2 emissions per capita. However, and similarly to the urban scaling case, the exponent $\alpha$ is also likely to be biased by the confounding effect of population, since density and population values are also correlated 24; 25.

In order to improve the description of urban emissions, we propose here an analogy with the economic theory of production functions. This theory has a central role in several branches of economics 26, and in general terms, a production function establishes a mathematical relationship between inputs (capital, labor, land, etc.) and output of goods (iron, cars, wheat, etc.) in some production process. Using this mathematical description, economists may ask how much output can be produced with particular combinations of inputs, and what are the alternatives (in terms of inputs) for producing a particular good. These ideas transpose well into our context assuming CO2 emissions as the output and population and area as the inputs in a production process mediated by cities. Similarly to a two-factor production model, we thus consider that $C=F(P,A)$, where $F(\dots)$ stands for the form of the production function. By putting this analogy forward, we establish a more general approach for modeling urban emissions that simultaneously accounts for the effects of population and area along with possible interactions between these urban metrics. As we shall verify, we can borrow not only functional forms from the theory of production functions but also key concepts that are very useful in the context of urban emissions (Methods).

Cobb-Douglas model of urban carbon emissions

We start this analogy with the Cobb-Douglas model 28, arguably the most widely known and used production function 26. In our case, it takes the form

[TABLE]

where $\beta_{P}$ and $\beta_{A}$ are two independent exponents. We immediately notice that this model recovers the urban scaling [Eq. (1)] if $\beta_{P}=\beta$ and $\beta_{A}=0$ (that is, when ignoring the effect of area) and the per capita density scaling [Eq. (2)] if $\beta_{P}=\alpha+1$ and $\beta_{A}=-\alpha$. We further remark that the Cobb-Douglas model can be obtained from Eqs. (1) and (2) if we consider the empirical relation between population and area 24; 25 (Methods). Similarly to the models of Eqs. (1) and (2), the Cobb-Douglas function exhibits a scale-invariant elasticity $\varepsilon=\beta_{P}+\beta_{A}$ (Methods), meaning that a proportionate increase in emissions associated with a proportionate increase in population and area is independent of $P$ and $A$. Thus, when $\beta_{P}+\beta_{A}<1$ there are decreasing returns to scale (doubling $P$ and $A$ implies less than doubling $C$), whereas if $\beta_{P}+\beta_{A}>1$ there are increasing returns to scale (doubling $P$ and $A$ implies more than doubling $C$), and only for $\beta_{P}+\beta_{A}=1$ this model presents constant returns to scale (doubling $P$ and $A$ implies exactly doubling $C$). Thus, we notice that Eq. (2) is a particular case of the Cobb-Douglas model with constant returns to scale. On the other hand, without any constraint for the exponents $\beta_{P}$ and $\beta_{A}$, the Cobb-Douglas model represents a genuine generalization that cannot be related to Eqs. (1) and (2). In addition to that, we can interpret Eq. (3) as the result of accounting for the confounding effect of area $A$ within the urban scaling framework [Eq. (1)] via a multiple linear regression (in log-transformed variables).

Although the model of Eq. (3) may represent a better description for CO2 emissions, it introduces some drawbacks related to the use of ordinary-least-squares for finding the best fitting parameters $\beta_{P}$ and $\beta_{A}$. This happens because population and area are correlated to each other, a problem known as multicollinearity and that can lead to unstable estimates for the model parameters. As detailed in Methods, we have applied the regularization approach of the ridge regression 29; 30 in order to account for this problem. To state briefly, the ridge regression adds a penalty/regularization term proportional to the square of the magnitude of coefficients upon the residual sum of squares, which in turn stabilizes the regression coefficients and accounts for the multicollinearity. This approach yields ${\beta}_{P}=0.31\pm 0.01$ and ${\beta}_{A}=0.45\pm 0.03$ ($p$-values $=0$, permutation test, Supplementary Figure 3), and Fig. 2A shows the relationship between the actual values of the CO2 emissions and those predicted by Eq. (3). We have verified that the Cobb-Douglas model provides a significantly better fit to our data (Supplementary Figure 2) when compared with the urban scaling [Eq. (1)] and the per capita density scaling approaches [Eq. (2)]. Moreover, the fact that ${\beta}_{P}$ is much smaller than $\beta$ reinforces the idea that the urban scaling approach is indeed affected by the confounding effect of the area.

Because $\beta_{P}+\beta_{A}<1$, our results indicate that CO2 emissions display diminishing returns as population and area are incrementally increased by the same factor (that is, keeping density constant). In particular, our estimates indicate that every $1$% increase in both the population and area of a city associates with a $0.76$% increase in its emissions. The interconnected role of population and area on CO2 emissions is better visualized in Fig. 2B, where we depict a contour plot of Eq. (3) on logarithmic scale. In this representation, the isoquants (or isolines) are described by straight lines [$\log A\sim-(\beta_{P}/\beta_{A})\log P$] and show how population and area must change to keep emissions constant. These isoquants also indicate that if the population of a city increases, its population density must also increase to keep emissions unaltered. This behavior is better understood by rewriting Eq. (3) as $C\sim P^{\beta_{P}+\beta_{A}}\,(P/A)^{-\beta_{A}}$ and noticing that $\beta_{P}>0$ and $\beta_{A}>0$ (Supplementary Figure 4 shows the contour plot in terms of population density). From this form of Eq. (3), we also conclude that a proportionate change in population has more impact on CO2 emissions than a proportionate change in density since $|\beta_{P}+\beta_{A}|>|\beta_{A}|$. For a particular urban unit in our data, this means that if its population decreases by 1% while its density remains constant (that is, a 1.01% increase in its area), the model predicts a 0.76% reduction in its CO2 emissions; whereas a 1% raise in its density while population is unaltered (that is, a 1% reduction in its area), implies only 0.45% decrease in its CO2 emissions.

Another interesting aspect of an isoquant is its slope, a quantity known as the technical rate of substitution in economics 26 (Methods). In our context, this slope measures how much the population of a city should change in response to alterations in its area in order to keep the same level of emissions. For Eq. (3), the slopes of the isoquants are $\frac{dA}{dP}=-\frac{\beta_{P}}{\beta_{A}(P/A)}$, and thus, they are completely determined by the city density (assuming that $\beta_{P}$ and $\beta_{A}$ are known). If we consider a logarithmic scale, these slopes are equal to $\frac{d\log A}{d\log P}=-\beta_{P}/\beta_{A}$ regardless of the values of $P$ and $A$ (as we see in Fig. 2B). In economics, these isoquants are also analyzed in terms of the so-called elasticity of substitution 26 (Methods), a dimensionless measure that (mapped to our case) quantifies the efficiency at which population and area substitute each other, and that somehow reflects the shape of the isoquants 26. Usually, more L-shaped isoquants are associated with low elasticity of substitution (that is, there is no room for replacing $A$ by $P$ while keeping emissions constant), whereas more linear/smooth isoquants tend to have high elasticity of substitution (it is easy to replace $A$ by $P$ while keeping emissions constant). The Cobb-Douglas model has unitary elasticity of substitution, regardless of the values of $P$, $A$, and $C$, and also the exponents $\beta_{P}$ and $\beta_{A}$ 26. Thus, although the Cobb-Douglas model provides a better fit to our data [compared with Eqs. (1) and (2)], it also makes a series of assumptions that do not have any compelling reasons to hold true in urban systems (as it also happens in economics 26). Moreover, Fig. 2A shows that Eq. (3) has a bias for large values of $C$, which indicates that relaxing some underlying assumptions of the Cobb-Douglas model may lead to a better description of the emissions.

We first relax the condition of unitary elasticity of substitution by considering a model based on the constant elasticity of substitution (CES) production function 31; 26

[TABLE]

where $\gamma$ is a parameter and $\beta_{P}+\beta_{A}=1$. The CES model emerged as a generalization of the Cobb-Douglas function (which is recovered when $\gamma\to 0$) exactly because economists have considered the assumption of unitary elasticity of substitution as unduly restrictive 31. As the name suggests, the CES model has a constant elasticity of substitution equal to $1/(\gamma+1)$, and thus, by varying $\gamma$ we have a wide range of possible elasticities. We have adjusted Eq. (4) to our data using the Levenberg-Marquardt algorithm. However, these fits are very problematic because of the non-linear nature of the model. Depending on the initial guess used for the parameters, we find very unstable and meaningless results. These fits also lead to very large variations in the parameters when applying a resampling strategy to our data. Furthermore, we have observed that even when the fits converge, the CES model does not represent the best description to our data. In addition to the non-linearity, this happens because the CES model also carries crucial assumptions that appear not to hold in our case. This model has unitary elasticity of scale ($\varepsilon=1$, that is, it assumes that by doubling $P$ and $A$ implies doubling $C$), and similarly to Cobb-Douglas, the slopes of the isoquants $\left(\frac{dA}{dP}=-\frac{\beta_{P}}{\beta_{A}(P/A)^{\gamma+1}}\right)$ are determined solely by the city density.

Transcendental logarithm model of emissions

To overcome these limitations and constraints, we considered a more general production function known as the transcendental logarithm (translog) model 32; 26. Just like the CES is a natural extension of the Cobb-Douglas function, the translog model represents the next logical step towards a more flexible function for modeling the CO2 emissions in terms of $P$ and $A$. This model is written as

[TABLE]

where $\beta_{C}$ is an additional parameter. Differently from Eqs. (3) and (4), the translog model has a non-constant elasticity of scale [$\varepsilon=\beta_{P}+\beta_{A}+\beta_{C}\log(PA)$], meaning that proportionate changes in the emissions associated with proportionate changes in population and area depend on the initial values of $P$ and $A$. Furthermore, the elasticity of substitution of Eq. (5) varies with $P$, $A$ and $C$, and the isoquant slope depends on $P$ and $A$ (Methods). By comparing the translog model with the Cobb-Douglas, we notice that all this additional flexibility is solely related to the inclusion of the interaction term [$\beta_{C}\log(P)\log(A)$] between population and area, where $\beta_{C}$ quantifies the intensity of this interaction. It is this interaction term that allows the effect of population and area on emissions to vary with $P$ and $A$. This effect is better understood in terms of the marginal products 26 (Methods), an economic quantity that mapped to our context represents the response in CO2 emissions caused by changes in population or area. The marginal product of population (in logarithmic scale) is defined as $\frac{d\log C}{d\log P}=\beta_{P}+\beta_{C}\log A$, while the marginal product of area is $\frac{d\log C}{d\log A}=\beta_{A}+\beta_{C}\log P$. Thus, for instance, we observe that the marginal product of population depends on the area of the urban unit. This behavior contrasts with the Cobb-Douglas predictions (corresponding to $\beta_{C}=0$), in which the marginal products are independent of $P$ and $A$.

As in the Cobb-Douglas case, we have considered the ridge regression approach in order to account for the multicollinearity and adjust Eq. (5) to our data (Methods). This approach yields $\beta_{P}=0.28\pm 0.02$, $\beta_{A}=0.14\pm 0.05$, and $\beta_{C}=0.07\pm 0.01$ ($p$-values $=0$, permutation test, Supplementary Figure 3). We further note that the translog model is a significantly better fit to our data when compared with the previously-discussed models (Supplementary Figure 2). This fact is highlighted in Fig. 2C, where we observe that the model of Eq. (5) refines the quality of the predictions and reduces the bias observed for Eq. (3) in urban areas with high emissions. We further remark that all fitted parameters of Eq. (5) are significantly different from zero, confirming the predictive power of the interaction term between population and area. Therefore, in addition to improving the description of the CO2 emissions, Eq. (5) reveals that the effect of population and area on the emissions intensifies with the increase of urban population and area. This becomes clear by noticing that the elasticity of scale increases with $P$ and $A$ [$\varepsilon=0.42+0.07\log(PA)$], that is, the more populous and the more widespread a city is, the larger is the impact of a proportionate change in its population and area (a growth with constant density) on its emissions. It is further intriguing to notice that, because the elasticity of scale varies with $P$ and $A$, the translog model displays decreasing, increasing, or constant returns to scale depending on whether the product $\Omega=PA$ is respectively smaller, larger or equal to the critical value

[TABLE]

For the US data $\Omega^{*}\approx 1.93\times 10^{8}$, and thus, cities having $\Omega<\Omega^{*}$ display decreasing returns to scale, while those having $\Omega>\Omega^{*}$ feature increasing returns to scale. For instance, the translog model predicts that a 1% increase in population and area of a large city with $P=8\times 10^{6}$ and $A=6,000$ km2 (roughly the size of Chicago, $\Omega=4.8\times 10^{10}$) associates with 1.17% raise in its emissions, whereas the same change in a relatively small city with $P=90,000$ and $A=140$ km2 (roughly the size of Santa Fe, $\Omega=12.6\times 10^{6}$) associates with only 0.92% raise in its emissions.

The interconnected role of population and area on CO2 emissions is better visualized in the contour plot of Eq. (5) shown in Fig. 2D. In comparison with Fig. 2B, we note that the interaction term bends the isoquants upward and make their slopes $\left[\frac{dA}{dP}=\frac{-1}{(P/A)}\left(\frac{\beta_{P}+\beta_{C}\log A}{\beta_{A}+\beta_{C}\log P}\right)\right]$ a function of $P$ and $A$, and not only of the city density as in the Cobb-Douglas model. We further observe that the spacing between the isoquants representing equally incremented values of $\log C$ changes with the values of $\log P$ and $\log A$. This behavior contrasts with the equally spaced isoquants produced by the model of Eq. (3), and emphasizes that proportionate changes in emissions caused by changes in both population and area depend not only on the intensity of the changes but also on the initial values of population and area.

In terms of density, we can rewrite Eq. (5) as

[TABLE]

and from this expression, we find that the elasticity of scale in terms of population and density is $\varepsilon=\beta_{P}+\beta_{C}\log A$. Thus, the more widespread a city is, the larger is the impact of a proportionate change of population and density (that is, area remains constant) on its emissions. In a concrete example, our estimates indicate that a $1\%$ raise in population and density of a large city with $6,000$ km2 associates with a $0.54\%$ increase in its emissions, while the same change in a city with $140$ km2 correlates with a $0.42\%$ raise in its emissions. Furthermore, in terms of population and density, increasing returns to scale ($\varepsilon>1$) is only possible for urban areas exceeding the critical value $A^{*}=10^{(1-\beta_{P})/\beta_{C}}$. For the US data, $A^{*}\approx 1.93\times 10^{10}$ km2, an area roughly corresponding to 38 times the area of Earth. Therefore, our estimates show that only decreasing returns to scale are possible when the CO2 emissions are described in terms of population and density.

The translog model of Eq. (7) also allows us to verify whether changes in population have more impact on the emissions than changes in density. To do so, we have compared the absolute values of the marginal products of population $\left[\frac{d\log C}{d\log P}=\beta_{P}+\beta_{A}+\beta_{C}\log(PA)\right]$ and density $\left[\frac{d\log C}{d\log(P/A)}=-\beta_{A}-\beta_{C}\log P\right]$. The marginal product of population represents the response in emissions associated with changes in population when density remains constant, whereas the marginal product of density expresses the response in emissions caused by a change in density when population remains constant. The inequality $\left|\frac{d\log C}{d\log P}\right|>\left|\frac{d\log C}{d\log(P/A)}\right|$ simplifies to $A>A^{*}=10^{-\beta_{P}/\beta_{C}}$ when $\beta_{P},\beta_{A},\beta_{C}>0$ and $A,P>1$ (assumptions that agree with our estimates). By plugging the estimated values of $\beta_{P}$ and $\beta_{C}$, this condition becomes $A>10^{-4}$ km2, and hence we conclude that population size always have more impact on emissions than changes in population density. We have reached the same conclusion with the simpler model of Eq. (3). However, and as we have verified, the translog approach further refines the description of CO2 emissions and indicates that the impact of population and density on the emissions changes with the population and density of the cities. For instance, according to our estimates for the US data, a 1% increase in the density of a town with 10,000 inhabitants associates with a 0.42% reduction in its emissions, while the same change in a larger city with 1 million inhabitants associates with a 0.56% reduction in its emissions. This behavior contrasts with results of Eq. (3), whose predictions related to changes in the density of a city are independent of its population. By carrying this predictions forward and in line with other studies 33, our results suggest that the densification of large populated urban areas is likely to have important contributions to the reduction of urban CO2 emissions.

Discussion

We have shown that two conventional approaches used to study the effect of urbanization on urban CO2 emissions suffer from confounding effects, and are unable to describe the interconnected role of population and area on urban emissions. Inspired by the economic theory of production functions, we have proposed new models for describing urban emissions simultaneously in terms of the population and area (or population density) of urban units. These models not only account for such confounding factors but significantly refine the description of the emissions in terms of urban quantities. In addition to being better fits to data, our models reveal intriguing aspects about the interplay between population and area (or density) on urban emissions that would be entirely neglected under the urban scaling or the per capita density scaling frameworks. Among these findings, our results indicate that variations in emissions associated with proportionate changes in population and area do not only depend on the magnitude of these changes but also on the product $\Omega=PA$. In particular, depending on whether $\Omega$ exceeds or not the critical value $\Omega^{*}$, urban emissions can display increasing ($\Omega>\Omega^{*}$) or decreasing returns to scale ($\Omega<\Omega^{*}$) with population and area. When described in terms of population and density, we have found that urban emissions display decreasing returns to scale, meaning that doubling population and density of a city always associates with less than doubling its emissions. We have further verified that changing the population of a city has more impact on its emissions than changing its density. In spite of that and in general terms, our models define conditions under which changes in population have more impact than changes in density (or vice versa) on emissions and further predict a transition-like behavior where the dominant role changes between these quantities if the urban area exceeds a threshold value ($A^{*}$).

Our work has however its limitations in a sense that ideally the comparison between the effects of population and area (or density) on the emissions should be made after accounting for every other factor (such as economic activity, technology, and even individual attitudes) that possibly affects urban emissions. Thus, while our models account for the confounding effects of area (and density), the emissions may also be affected by other confounding factors not available in our dataset. One possibility for addressing this problem would be to include further control variables in our models, an approach that somehow resembles the IPAT equations 34; 35; 36, a framework proposed to model environmental impact (I) as the product of population (P), affluence (A), and technology (T), but with the advantage of considering population density (or area) as a predictor and allowing the interactions among such factors. Another possibility for overcoming these possible confounding effects is to combine our approach with the recently proposed urban Kaya scaling 20 that relates CO2 emissions, population, gross domestic product, and energy consumption. Combining these different approaches into a single and coherent framework could represent an exciting perspective for solving the economics of urban CO2 emissions and defining its most important covariates. However, such endeavors require homogeneous and consistent data, which are still scarce on large spatial scales. While moving from urban units defined in terms of connected urban spaces to some political or administrative divisions would be a possibility, this approach is likely to introduce serious bias to the empirical estimates 37 in addition to overestimating urban areas 38 (Methods). Another important limitation of our study is related to the intra-city processes and urban characteristics that cannot be accounted for only by population and area (or density). Case studies on this subject have shown that the urban form and intra-city population distribution have a substantial impact on urban emissions, particularly on transportation emissions. Cities from rapid developing countries such as China and India have undergone through a remarkable decentralization and suburban growth processes 39. These more dispersed urban forms and the consequent increase of the population living in urban frontier areas have direct implications for commuting and contribute to increasing CO2 emissions 40; 39; 41. Regarding this aspect, it would be very interesting for future works to include possible covariates able to account for population imbalance and urban form in our models and thus quantify their impact in a large scale study.

Despite of these limitations, our work adds to the current understanding about the role of urbanization on CO2 emissions, shedding light particularly upon the role of population and urban area (including their interactions) on urban emissions. Such interactions are completely overlooked within the urban scaling and per capita density scaling approaches and our work demonstrates that they play an important role on the description of urban emissions. Finally, our framework can be directly applied to other urban metrics in the place of emissions, opening thus a considerable range of possibilities for investigating the interplay between population and area (or population and density) over other important urban metrics.

Methods

Dataset

Our dataset is the same as analyzed by Gudipudi et al. 15 and comprises the CO2 emissions from urban areas of the United States (US) in the year 2000. As described by Gudipudi et al., this dataset is compiled from different sources through the following steps. Firstly, gridded population data are obtained from the Global Rural-Urban Mapping Project (GRUMP) 42 and the Global Land Cover Dataset (GLC) 43. Both datasets are from the year 2000 and are available at a spatial grid resolution of 1 km $\times$ 1 km. The GRUMP data are spatially overlaid to the GLC in order to attribute population to the land use, which is classified as urban and non-urban. Next, sectoral emissions data (building and transportation) are obtained from the Vulcan Project 44. This dataset is initially available at a resolution of 10 km $\times$ 10 km and has been down-scaled to 1 km $\times$ 1 km to be superimposed on the populated settlements. This process consists of equally splitting the emissions located in a cell of the Vulcan project among all overlapping population cells classified as urban. Finally, the city clustering algorithm (CCA) 27 is used to systematically define the urban units, leading to the population size ($P$ in raw counts), area ($A$ in square kilometers), and CO2 emissions ($C$ in tonnes of CO2) for each urban unit. CCA is an iterative clustering algorithm that assigns any two cells to the same cluster if their distance is smaller or equal than a predefined threshold distance $l$.

All results presented in our manuscript have been obtained using $l=5$ km, a threshold distance that does not overestimate nor underestimate the urban extents 15. However, we have verified that our conclusions are very robust against variations in $l$ from $l=1$ km to $l=10$ km. In particular, the translog model [Eq. (5)] always provides the best description for the US data (Supplementary Figure 5). We have observed that the scaling exponents $\beta$ and $\alpha$ present slightly variations with threshold distance $l$ (Supplementary Figure 6) that have no implications for the conclusions drawn from our findings. The parameters of the Cobb-Douglas ($\beta_{P}$ and $\beta_{A}$) and translog ($\beta_{P}$, $\beta_{A}$, and $\beta_{C}$) moldels present somewhat larger variations (Supplementary Figure 7), but all remain statistically significant, particularly the parameter related to the interaction term in the translog model ($\beta_{C}$). These changes affect our point estimates for the Cobb-Douglas (Supplementary Figure 8) and translog (Supplementary Figure 9) models. In general, changes in emissions associated with a 1% increase in $P$ and $A$ tend to increase with $l$. Similarly, the reduction in emissions associated with -1% change in $P$ with fixed density or with a 1% change in density with fixed $P$ also increase with $l$. This dependence on $l$ is smaller in the translog than in the Cobb-Douglas model. The critical product $\Omega^{*}$ displays a decreasing trend with $l$ (Supplementary Figure 10), that in turn affects the point from which the effect of $P$ and $A$ on $C$ changes from decreasing to increasing returns to scale. The critical value $A^{*}$ is approximately independent of $l$ (Supplementary Figure 10), indicating that changes in population have more impact than changes in density on emissions regardless the value of $l$. We have further verified that the contour plots of the translog function [Eq. (5)] are very similar for different values of $l$ (Supplementary Figure 11).

These robustness tests are important because there has been a great debate about how to accurately define the correct boundaries of a city 37. In spite of that, there is still no consensus on this issue nor has a fail-safe procedure for defining the correct boundaries of a city been proposed yet. This issue also has great similarity with the more general concept of clustering and a quite similar issue arises when applying community detection algorithms in complex networks. All these topics have been exhaustively studied, but no silver bullet method exists. In the case of cities, additional complexity emerges because some urban indicators are more spatially constrained than others, and also because people commute to work and move from place to place in the long run.

Partly because of the seminal works by Bettencourt et al. 21, the use of functional definitions for cities such as the Metropolitan Statistical Areas (MSAs) in US or larger urban units (LUZ) and metropolitan areas (MAs) in Europe have become very popular. These definitions are based on the idea of integrated socio-economic units and appear to be the gold standard for the urban scaling hypothesis as well as other purposes. In particular, MSAs are defined by a core county (or even more than one) having at least 50,000 people aggregated with adjacent counties that display a high degree of interaction (social and economic) with the central county (as measured by commuting flows). While this definition may work well for studying urban scaling, it is very problematic in our case. Since MSAs are made up of counties, they often include vast rural areas which in turn hugely overestimate the urban extent areas. In addition to that, MSAs can also fragment urban clusters into different pieces.

These problems are the main reason why we have chosen the CCA to define the urban units in our study. It is worth noticing that the CO2 emissions we have analyzed are from building and transportation, and thus primarily associated with settlements where people reside and commute. Because of that, we argue that the approach of combining gridded data from the Vulcan project, GLC and GRUMP with CCA (proposed in Ref. 15) provides more precise emission estimates and associated spatial extents of urban clusters. Despite these problems and limitations, we have also applied our models to emissions data associated with MSAs. To do so, we have used the dataset provided by Fragkias et al. 8 and considered the emissions from the year 2000 (the same used in our analysis). Indeed, these data allowed us to consider not only MSAs but also Micropolitan Areas ($\mu$SAs are defined as labor market areas with population between 10,000 and 50,000 people that are also made up of counties) and both together (Core Based Statistical Areas – CBSAs).

Supplementary Figure 12 shows the urban and per capita density scaling laws for MSAs. We notice that the quality of these relationships is not comparable with those reported in Fig. 1, and a similar situation happens for $\mu$SAs and CBSAs. Supplementary Figure 13A compares the scaling exponents obtained from MSAs, $\mu$SAs, and CBSAs with those obtained via CCA with $l=5$ km. We observe that the values of $\beta$ estimated from these functional city definitions are much closer to one (in agreement with Fragkias et al. 8) than the values obtained with the CCA approach. As discussed in more detail by Bettencourt et al. 45, the disaggregation from the true urban unit can either introduce a positive or negative bias in the estimates of $\beta$. On the other hand, the aggregation of different true urban units tends to make $\beta$ closer to one. In the case of MSAs, it is likely that both disaggregation and aggregation effects play some role, but the fact that $\beta$ is smaller for $\mu$SAs than MSAs suggests that aggregation may have greater influence. Due to the poor quality of the per capita density scaling, it is hard to directly compare the values of $\alpha$ obtained from MSAs, $\mu$SAs and CBSAs with those estimated with the CCA; however, Supplementary Figure 13B shows that at least they have the same sign.

We have also applied the models of Eqs. (3) and (5) to MSAs, $\mu$SAs, and CBSAs data. Supplementary Figure 14 shows that these models do not represent an improved description when compared with Eq. (2). For CBSAs and MSAs, Eq. (3) has about the same predictive power as Eq. (2). Supplementary Figure 15 compares the exponents $\beta_{P}$ and $\beta_{A}$ obtained from MSAs, $\mu$SAs and CBSAs with the CCA values for $l=5$ km. We notice that $\beta_{P}$ is not so different from $\beta$ for these functional city definitions, indicating that the confounding effect of the area is much weaker when compared with the CCA results. This fact reinforces the idea that MSAs and $\mu$SAs areas are not good predictors for CO2 emissions. We also observe that $\beta_{P}$ is much larger for the functional definitions, and that $\beta_{A}$ is smaller than the values obtained from the CCA approach. In spite of all discrepancies, the results for MSAs, $\mu$SAs and CBSAs also indicate that population has more impact on emissions than density because $|\beta_{P}+\beta_{A}|>|\beta_{A}|$. On the other hand, $\beta_{P}+\beta_{A}\approx 1$ for the functional definitions, while $\beta_{P}+\beta_{A}<1$ for the CCA approach. The approximate constant returns to scale observed for the functional city definitions is also likely to be related to disaggregation and aggregation effects that we previously discussed.

Finally, it is worth remarking that the CCA is also likely to suffer from disaggregation or aggregation effects, as is the case of administrative or functional city definitions. However, differently from such ad hoc definitions, CCA allows us to quantify the impact of such effects by changing the threshold distance $l$ and to verify that our conclusions are robust under different values of $l$.

Cobb-Douglas and the urban scaling models

To relate the Cobb-Douglas model [Eq. (3)] with the urban scaling [Eq. (1)] and the per capita density scaling [Eq. (2)], we first rewrite Eq. (1) as

[TABLE]

where $\beta=\beta_{1}+\beta_{2}$ stands for the same urban scaling exponent as in Eq. (1). Now, we replace the right-most $P$ in the previous equation by the allometry relation between population and area 24; 25, $P\sim A^{\delta}$ (where $\delta$ is another power-law exponent), leading to

[TABLE]

which has the same form as Eq. (3), that is, $\beta_{1}=\beta_{P}$ and $\beta_{2}\delta=\beta_{A}$. Consequently,

[TABLE]

and the two approaches are uniquely related only if there is another constraint for the exponents. One possibility is by imposing constant returns to scale in the Cobb-Douglas model ($\beta_{A}+\beta_{P}=1$), which leads to

[TABLE]

for $\delta\neq 1$; while for $\delta=1$, we obtain $\beta_{P}=1$ and $\beta_{A}=0$. Thus, the existence of an additional constraint for the parameters $\beta_{P}$ and $\beta_{A}$ implies that the Cobb-Douglas model is equivalent to Eqs. (1) and (2); otherwise, it represents a generalization.

The Cobb-Douglas model with $\beta_{A}+\beta_{P}=1$ is also related to scaling relationships between indicator density and population density 46; 47; 48. To obtain this connection, we rewrite Eq. (3) as

[TABLE]

where $\beta_{P}=\theta$ and $\beta_{A}=1-\theta$. Next, we divide both sides by $A$

[TABLE]

leading to a scaling relationship between CO2 density and population density.

Analogy with the production functions

As we have argued, our approach is inspired by the economic theory of production functions 26. By following this analogy, we have considered the urban emissions as the output and population and area (or density) as the inputs of a two-factor production process mediated by cities. The mathematical formula that describes the possible relations between the inputs ($P$ and $A$) and the output ($C$) is the production function, that is, $C=F(P,A)$. The functional forms for $F$ used in our work comprise the most widely known and used production functions in economics and should be viewed as an empirical/phenomenological description (as it also happens in economics). In what follows, we summarize concepts from the economic theory of production functions that have been used in our work.

Elasticity of scale. The elasticity of scale $\varepsilon$ is the ratio between a proportionate change in the output (emissions) and a proportionate change in the inputs (population and area), that is, $\varepsilon=(dF/F)/(dP/P)$, where $dP/P=dA/A$ represents the proportionate change in the inputs. This measure quantifies the impact of changing population and area on emissions.

Technical rate of substitution. The technical rate of substitution measures the rate at which an input must change in response to a change in the other input so that the output remains constant. In absolute value, it represents the slope of the isoquants of the production function. For instance, assuming a particular value for the output $F(P,A)=C_{0}$, the technical rate of substitution between area and population is $\delta_{P}^{A}=dA/dP$.

Elasticity of substitution. The elasticity of substitution $\sigma$ somehow summarizes the shape of an isoquant. It is defined as the ratio between a proportionate change in the inputs and the associated proportionate change in the slope of the isoquant. Mathematically, we write $\sigma=\frac{[d(P/A)/(P/A)]}{[d(dP/dA)/(dP/dA)]}=\frac{d\log(P/A)}{d\log(dP/dA)}$, where the numerator represents the proportionate change in the inputs and the denominator the proportionate change in the slope of the isoquant. This measure quantifies the efficiency at which population and area substitute each other.

Marginal products. The marginal product of an input is defined as the (infinitesimal) change in the output resulting from a (infinitesimal) change in one of the inputs. For instance, the marginal product of population is $\Delta_{P}=dF/dP$ (it can also be defined in terms of logarithmic quantities: $\Delta_{P}=d\log F/d\log P$).

We summarize all these properties calculated for Cobb-Douglas [Eq. (3)] and translog [Eq. (5)] models in Supplementary Table 2.

Fitting models with the ridge regression approach

As we have discussed in the main text, multicollinearity is present in the models of Eqs. (3) and (5). This effect happens when at least two predictors in a multiple linear regression are correlated to each other 30; 29. Under this situation and depending on the degree of correlation among the predictors, ordinary least-squares estimates of the parameters can be unstable against minor changes in the input data and also display large standard errors. To better illustrate this problem, consider the simple linear model

[TABLE]

where $y$ is the response variable, $x_{1}$ and $x_{2}$ are the predictors, and $a_{1}$ and $a_{2}$ are the linear coefficients. The least-squares estimator for the parameters is usually written as $\bm{a}=\left[\begin{smallmatrix}a_{1}\\ a_{2}\end{smallmatrix}\right]=(\bm{X}^{T}\bm{X})^{-1}\bm{X}^{T}\bm{y}$, where $\bm{y}=\left[\begin{smallmatrix}y_{1}^{(1)}\\ \vdots\\ y^{(n)}\end{smallmatrix}\right]$ is a $n\times 1$ vector of the response variables, $\bm{X}=\left[\begin{smallmatrix}x_{1}^{(1)}~{}x_{2}^{(1)}\\ \vdots\hskip 17.07182pt\vdots\\ x_{1}^{(n)}~{}x_{2}^{(n)}\end{smallmatrix}\right]$ is a $n\times 2$ matrix of the regressors, and $n$ is the number of observations. If the values of predictors are strongly correlated, the inversion of the matrix $\bm{X}^{T}\bm{X}$ can become unstable, and consequently lead to unstable estimates for the linear coefficients.

To account for the multicollinearity problem, we have fitted Eqs. (3) and (5) by using the ridge regression approach 30; 29. This method solves the matrix inversion problem by adding a constant $\lambda$ to the diagonal elements of $\bm{X}^{T}\bm{X}$, so that the ridge estimator for the linear coefficients is $\bm{a}=(\bm{X}^{T}\bm{X}+\lambda\bm{I})^{-1}\bm{X}^{T}\bm{y}$, where $\bm{I}$ is the identity matrix. The ridge estimation is equivalent to finding the optimal linear coefficients that minimize the residual sum of squares plus a penalty term (also called regularization parameter) proportional to the sum of the squares of the linear coefficients 30; 29, that is, finding the $\bm{a}$ that minimizes the objective function $\lVert\bm{y}-\bm{X}\bm{a}\rVert^{2}+\lambda\lVert\bm{a}\rVert^{2}$. The optimal value of $\lambda$ is usually unknown in practice and needs to be estimated from data. To do so, we have used the approach of searching for the value of $\lambda$ that minimizes the mean squared error (MSE) in a leave-one-out cross validation strategy. In this approach, we estimate $\bm{a}$ (for a given $\lambda$) using all data except for one point that is used for calculating the squared error. This process is repeated until every data point is used exactly once for estimating the squared error, and then we calculate the value of the MSE for a given $\lambda$. The optimal value of $\lambda=\lambda^{*}$ is the one that minimizes the average value of the MSE estimated with the leave-one-out cross validation method. We have also standardized all predictors before searching for the optimal value $\lambda^{*}$. This is a common practice when dealing with regularization methods and ensures that the penalty term is uniformly applied to the predictors, that is, the normalization makes the scale of the predictors comparable and prevents variables with distinct ranges from having uneven penalization.

The standardized version of Eq. (3) can be written as

[TABLE]

where

[TABLE]

In addition, $\mu_{P}$ is the average value of $\log P$ ($\mu_{P}=\frac{1}{n}\sum\log P$), $\sigma_{P}$ is the standard deviation of $\log P$ [$\sigma^{2}_{P}=\frac{1}{n-1}\sum(\log P-\mu_{P})^{2}$], $\mu_{A}$ is the average value of $\log A$ ($\mu_{A}=\frac{1}{n}\sum\log A$), and $\sigma_{A}$ is the standard deviation of $\log A$ [$\sigma^{2}_{A}=\frac{1}{n-1}\sum(\log A-\mu_{A})^{2}$]. It is worth remarking that the values of $\widehat{\log P}$ and $\widehat{\log A}$ are invariant against changes in the scale of $P$ and $A$, that is, their values do not change under the transformations $P\to\nu P$ and $A\to\nu A$, where $\nu$ is a positive constant. The same invariance holds for $\sigma_{P}$ and $\sigma_{A}$, whereas the average values change according to $\mu_{P}\to\log\nu+\mu_{P}$ and $\mu_{A}\to\log\nu+\mu_{A}$.

The connection between the parameters of the standardized model ($\tilde{\beta}_{P}$ and $\tilde{\beta}_{A}$) and the usual ones ($\beta_{P}$ and $\beta_{A}$) is obtained by plugging Eq. (16) into Eq. (15), collecting the terms multiplying $\log P$ and $\log A$, and then directly comparing the results with Eq. (3). By following this approach, we find that

[TABLE]

where we observe that $\beta_{P}$ and $\beta_{A}$ are independent of the units of $P$ and $A$ (as they should be, since $\sigma_{P}$ and $\sigma_{A}$ are scale invariants). Supplementary Figure 16 illustrates how the ridge regression approach is applied to the model of Eq. (15). Supplementary Figure 16A shows the dependence of the MSE on the values of $\lambda$, from which we obtain $\lambda^{*}=9.78$. Supplementary Figure 16B shows the dependence of the parameters $\tilde{\beta}_{A}$ and $\tilde{\beta}_{P}$, whose values for the optimal regularization term ($\lambda=\lambda^{*}$) are $\tilde{\beta}_{P}=0.37\pm 0.02$ and $\tilde{\beta}_{A}=0.24\pm 0.02$, which correspond to ${\beta}_{P}=0.31\pm 0.01$ and ${\beta}_{A}=0.45\pm 0.03$. The errors in $\tilde{\beta}_{P}$ and $\tilde{\beta}_{A}$ stand for the standard deviation of their values estimated over 1000 random samples with replacement of the data, as shown in Supplementary Figure 17. The errors in ${\beta}_{P}$ and ${\beta}_{A}$ are calculated with common error propagation formulas. Moreover, the $p$-values of permutation tests reject the null hypothesis that these parameters are equal to zero.

In the case of Eq. (5), its standardized version can be rewritten as

[TABLE]

where the connecting formulas

[TABLE]

are obtained as in the previous case. Unlike the models of Eqs. (1), (2), and (3), the generalization expressed by Eq. (5) is not scale-invariant and thus its parameters depend on the measurement units. In particular, if the area is rescaled by a factor $\nu$ ($A\to\nu A$), $\beta_{P}$ is incremented by the factor $-\beta_{C}\log\nu$ and $\beta_{A}$ remains unchanged. Similarly, if population is rescaled by $\nu$ ($P\to\nu P$), $\beta_{A}$ is also incremented by $-\beta_{C}\log\nu$ and $\beta_{P}$ remains unchanged. Only $\beta_{C}$ is invariant against scale changes in $P$ and $A$. Thus, all interpretations related to the behavior of the emissions obtained from the model Eq. (5) involve the assumption that area is expressed in units of km2 and raw population counts. Supplementary Figure 18 illustrates how the ridge regression approach is applied to the model of Eq. (18). Supplementary Figure 18A shows the dependence of the MSE on $\lambda$, from which $\lambda^{*}=8.67$ is obtained. Supplementary Figure 18B shows the behavior of the model parameters as a function of $\lambda$. This approach yields $\tilde{\beta}_{P}=0.40\pm 0.02$, $\tilde{\beta}_{A}=0.17\pm 0.02$, and $\tilde{\beta}_{C}=0.044\pm 0.006$ for $\lambda=\lambda^{*}$, which correspond to $\beta_{P}=0.28\pm 0.02$, $\beta_{A}=0.14\pm 0.05$, and $\beta_{C}=0.07\pm 0.01$. The standard errors are calculated as in the previous case and the $p$-values of the permutation tests reject the null hypothesis that the model parameters are equal to zero (Supplementary Figure 3).

In addition to the models of Eqs. (3) and (5), we have further tested for a more general translog model (full translog model) having quadratic terms in $\log P$ and $\log A$, that is,

[TABLE]

where $\beta_{P}^{\prime}$ and $\beta_{A}^{\prime}$ are additional parameters. This expression can also be related to the CES production by applying the Taylor series expansion to Eq. (4) around the point $\gamma=0$ (the Kmenta approximation 49). We have fitted Eq. (20) by following the same procedure used for Eqs. (3) and (5). In particular, Supplementary Figure 19 illustrates how the ridge approach is applied to this model and Supplementary Figure 20 shows the best fitting parameters and their estimated errors. However, as shown in Supplementary Figure 2, the full translog model of Eq. (20) does not improve the goodness of the fit when compared with Eq. (5).

Data availability

The dataset used in this study were obtained from Gudipudi et al. 15, which in turn rely on freely available data obtained from the Global Rural-Urban Mapping Project (GRUMP) 42, the Global Land Cover Dataset (GLC) 43, and the Vulcan Project 44. All data supporting the findings of this study are available from the corresponding authors on reasonable request.

Code availability

The code used for analyzing data is available from the corresponding authors on reasonable request.

Author contributions

H.V.R, D.R., J.P.K. designed research, performed research, analyzed data, and wrote the paper.

Competing interests

The authors declare no competing interests.

Acknowledgements.

We thank R. Gudipudi for providing the data and useful discussions. This work emerged from ideas discussed at the symposium Cities as Complex Systems (Hanover, July 13th-15th, 2016) which was generously funded by Volkswagen Foundation. H.V. Ribeiro thanks the financial support of CNPq (grants 303250/2015-1, 407690/2018-2, and 303121/2018-1). D. Rybski is grateful to the Leibniz Association (project IMPETUS) for financially supporting this research. J.P. Kropp acknowledges financial cooperation support (Brazil: KoUP1 – Remus) by University of Potsdam.