Climate is commonly defined as the long-term weather conditions over a specific region. These longterm conditions are usually approximated by the statistical description in terms of central tendency and variability of weather over a long period of time (i.e., 30 years) 1 . In particular, the average of weather variables over such a period is referred to as climatology or climate normal and, in practice, it is frequently used as a shorthand description of climate. However, this description is incomplete as it ignores variability and omits information that may be relevant for a variety of users/applications. The same is true for using only changes in the mean to describe climate change and to estimate its impacts, which are known to be a function of both changes in central tendency measures and variability [2][3][4] . Integrated assessment models of the climate and the economy usually translate changes in global or regional mean annual surface temperatures to economic damages through stylized mathematical representations, called damage functions 5,6 . Temperature change in IAMs is assumed to be gradual and smooth, and excludes any variability in the time and spatial domains 7,8 . This omission has two important implications. First, intra-and interannual variability are assumed to have no effect on the warming rate. All months and seasons warm by the same amount and variability modes, such as the Atlantic Multidecadal Oscillation (AMO) and the Pacific Decadal Oscillation (PDO), cannot amplify nor dampen the warming trend, nor is there any interplay between variability and climate change that affect the projected damages. These are oversimplifications that are useful for keeping IAMs as simple and transparent as possible, but do not fare well with what is known about climate change: warming rates do vary between months and seasons 9,10 ; low-frequency oscillations modulate the warming rate; there are important interactions and synergies between natural variability and the response of the climate system to changes in external forcing, and natural variability affects the magnitude of damages 8,11 . Second, in IAMs warming is incorrectly assumed to be spatially homogeneous over the globe and/or world regions. Observations, as well as climate models' projections, show that warming rates and changes in other climate variables (e.g., precipitation) vary widely over space 12,13 . In IAMs, global (or, in some cases, regional) annual temperatures act as an index of changes in the relevant climate variables over space (and time) which drive the physical and economic losses that damage functions aim to represent. While there is considerable evidence in favor of spatial stationarity of changes in climate, this does not imply that the spatial variance of changes remains constant for different levels of warming 2,14 . On the contrary, differential warming around the globe is clear from all global circulation climate models and implies that the variance of changes in temperature increases with global temperature rise. As an illustration, consider the case of the pattern scaling technique which has been extensively evaluated and is consider to adequately approximate the output from complex physical climate models, particularly in the case of mean surface temperatures [14][15][16] . This technique assumes spatial stationarity of changes in variables such as temperatures and precipitation. Under such assumption, a matrix of fixed weights 𝑃 𝑣 = 𝑝 𝑖,𝑗 𝑣 , which represents the spatial pattern of change for variable 𝑣 and for each pair of geographical coordinates 𝑖, 𝑗, can be obtained from climate models' output to approximate what the magnitude of change at the grid cell would be for variable 𝑣 and for a given value of change in global annual temperature ∆𝑇 𝑔 . Even in this case of spatial stationarity, the variance of the pattern 𝑃 𝑣 would increase as warming rises due to the properties of the variance: 𝑣𝑎𝑟(𝑣) = 𝑣𝑎𝑟(∆𝑇 𝑔 𝑃 𝑣 ) = (∆𝑇 𝑔 ) 2 𝑣𝑎𝑟(𝑃 𝑣 ). As is shown in the next sections, similar results are obtained using warming output from global physical models. Having spatial climate change projections that are consistent with physics is important for generating economic impacts of climate change in recently developed local IAMS in which the economic exposure to climate impacts and vulnerability varies across space 17 . Few studies have focused on incorporating the effects of climate variability or stochasticity into commonly used climate damage functions. These modified functions aim to represent the impacts of climate change in terms of mean and variability, which includes intra-annual to decadal, centennial, and even longer frequency oscillations. However, these few studies have been restricted to include variability in the time domain, ignore spatial variation, and do not account for the seasonality of climate impacts. Here we show that these limitations can have significant effects on their estimates and conclusions. Estrada and Tol 7 proposed recalibrating existing damage functions in such a way that their parameters account for the stochasticity in climate time series projections. They focus on global average temperatures and conclude that the mean monetary impacts produced by the stochastic damage function is similar those of the deterministic function, but that the welfare impacts are larger 7 . A recent study 18 suggests that accounting for interannual variability in global mean temperature time series, the risk of economic damages increases by trillions of dollars even if variability is assumed to be moderate. However, the expected losses are not biased and are the same as when temporal variability is not included. They use a zero-dimensional stochastic climate model and the Weitzman damage function, one of the damage functions with higher response to temperature change, to project the distribution of losses 18 . The authors calculate the risk premium of the effects of including internal variability in the time domain and combined with the uncertainty of the equilibrium climate sensitivity. They find that while there is no considerable effect on the estimated expected losses, the risk premium can be up to $46 trillion dollars (or about 58% of global output).
Of this risk premium about US$32 trillion areis due to aleatory uncertainty and US$14 trillion to the interaction of aleatory and epistemic uncertainty. These studies provide relevant information about the effects of internal variability in the time domain on economic impacts and conclude that omitting this factor can have important effects on some aspects of the estimation of the consequences of climate change. Another study looking at past impacts showed that there are interaction effects between natural variability and external forcing that could bias the estimated costs of climate change 8 . That is, variability may affect not only risk premium but also the estimated level of economic losses. In the present study, we show analytically and numerically that the way internal variability has been handled in current studies of the economic costs of climate change leads to important downward biases. We expand previous studies by considering temporal inter-and intra-annual variability and spatial variation, as well as by developing new annual, seasonal and monthly damage functions designed to include climate variability and take advantage of sectoral estimates of damages available in the literature. Our results show that considering internal climate variability, both in the time and space domains, significantly increases the projected economic damages, not only the risk premium a social planner would pay to avoid stochasticity. Therefore, the current estimates importantly underestimate the economic costs of climate change. Furthermore, we provide an analytic solution to include climate variability that does not require conducting simulation experiments and that is embedded in the most common currently used damage functions. This extension can be easily adopted in most IAMs currently available, as we illustrate by extending the RICE damage functions 19,20 to calculate the regional social cost of carbon. To include temperature change variability in time and space we incorporate three extensions to the commonly used climate-economy IAM framework (see Methods). First, coupled ocean-atmosphere general circulation models' output is used to create the monthly and spatially explicit temperature scenarios. Two approaches are used for this extension, one based on monthly scaling patterns (SP) 14,16,21 calculated by means of linear regressions (Methods) to produce gridded temperature projections from the output of MAGICC7 22 , a reduced complexity climate model. The second Figure 1 compares the economic impacts with and without climate variability. Moreover, it decomposes these results into spatial and temporal effects. This figure underlines the important role of Svar in the change in temperature for projecting the economic damages of climate change. Allowing for both spatial variation and temporal variability (S&Tvar) in ∆𝑇 𝑡 leads at the end of this century to damages of about 9.65% of global GDP, which are about 2.0% of global GDP higher than when no variability is included in the calculations (7.66% of global GDP). The largest change in the projected damages can be attributed to Svar, with damages raising to 9.11% of global GDP. The effect of incorporating Tvar by itself is quite small (0.01% of global GDP) and similar to what has been reported previously in the literature 7 . When variability in both spatial and temporal domains is allowed, the interaction effects of these two domains makes losses noticeably larger, with an increase in damages of about 0.54% of global GDP. The differences in damages increase with warming: according to the SP estimates, by 2050 if S&Tvar is ignored, damages would be underestimated by about US$1.6 trillion (US$1.4-US$2.3 trillion; Table S3) and about US$20 trillion by 2100 (US$17-US$28 trillion; Table S3). These results are like those obtained with the RM approach, except for those obtained with the CanESM5 which projects about 2ºC larger than the CMIP6 ensemble mean by 2085. Also, since the RM estimates are based on 30-year rolling statistics, the last year that can be calculated is 2085 and not 2100. In consequence, the magnitudes of these estimates are typically lower than those of the SP for 2100. Figure 2 shows the present value of the losses accumulated during the period 2020-2100 for each of the annual damage functions, using a 4% discount rate. Tables S3-S4 presents a sensitivity analysis using other climate models and both the SP and RM approaches. The present value of losses obtained with the damage function that excludes temperature change variability amounts to US183,841 billion, while the present value including S&Tvar is US$231,363 billion (US$225,679-US$249,750; Table S4), which represents a difference of about 55% of global GDP in 2020 (US$47 trillion with a range of US$42-US$66 trillion; about 26% higher than the original damage function). The present value of losses for the Svar damage function is US$218,489 billion (US$210,443-US$224,870), and US$184,086 (US$183,880-US$184,683) for Tvar which is much closer to that of the original damage function. The results in Table S3 show that not only the differences in the level of warming by various models (RM) strongly influence the projected losses, but that even when the differences in warming are controlled for (i.e., using the pattern scaling method; SP), the variability in space and time of each model increases damages between 23%-36%. Decomposing these additional losses (Table S4) reveals that about 65% (62%-75%) of them can be attributed to spatial variation and about 1% (-3%-2%) to temporal variability. The interaction between spatial variation and temporal variability has a large effect, accounting for about 33% (24%-45%) of the increases in losses. The effects of incorporating variability in the space and time domains are also important when comparing the costs associated to different reference and policy emissions trajectories. As shown in Table S5, when S&Tvar is included, the estimated damages are considerably larger for any emissions scenario. Moreover, using the SSP585 as an illustration, Table S5b-c shows that the benefits (in terms of avoided damages) are also larger for all levels of mitigation effort deemed possible under the SSP5 socioeconomic scenario 31 . When using the monthly damage function with all sectors aggregated by month, the calculated damages are very similar to those obtained with the annual damage function (Tables S3-S5).
For instance, by 2100 using the S&Tvar damage function, the present value of losses is about 0.03% of global GDP smaller compared with that of the annual damage function. These differences are due to the heterogeneity in warming during the year. The differences become clearer when monthly damages are analyzed per sector (Tables S6-S11). While including spatial variation always leads to higher damages, temporal variability can lead to lower damages in some cases when compared with projections from the annual damage function. These are sectors that are strongly influenced by seasonal variation, such as agriculture, forestry, water, and air pollution (Figure S3, S4). This occurs because the average annual warming is higher than the increase in temperatures in the months that are most relevant for the sector. In contrast, for some sectors temporal variability increases impacts significantly, as the warming that occurs in the relevant months is larger than the annual mean. These sectors are amenity, energy, health and time use, a sector for which warming and seasonal variability produces benefits. This illustrates the importance of relaxing the assumption of zero variability that is implicitly made in most of the current damage functions (equation 1). As suggested by Figure S4, damages on health, amenity and water would occur earlier in this century, reaching monthly losses of exceeding 0.1% of global GDP during the 2050s and doubling in the 2070s. If mitigation efforts that lead to an emissions trajectory like that of the SSP245 are achieved, exceeding 0.1% of GDP in monthly damages in sectors such as health could be delayed to the end of the 2080s, and losses larger than 0.2% would not be attained during this century. Figure 3 shows the present values calculated per sector under the SSP5-8.5, using monthly damage functions and different assumptions about temperature change variability (see Tables S6-S11 for other SSP scenarios and climate models). Except for Time use, climate change has negative effects in all sectors over this century and their magnitudes are highly dependent on the assumptions made about the temporal variability and spatial variation of warming. The largest losses correspond to the S&Tvar damage functions. Compared with the losses obtained under the no-variability assumption, the aggregated present value is about US$47 trillion larger and the sectors with larger differences are Health (US$12 trillion), Amenity (US$10.5 trillion), Catastrophe (US$6.3 trillion) and Time use for which climate change represents considerable benefits (US$9.1 trillion). In the case of the Svar damage function, the aggregated present value is about $35 trillion larger when compared with the damage function that assumes no variability. Although the magnitudes are smaller, the sectors with larger losses are the same: Health ($6.6 trillion), Amenity ($5.8 trillion), Catastrophe ($4.6 trillion) and Time use (-$5.0 trillion). The differences in present values between Tvar and the original novariability damage functions are much smaller, about $813 billion. Tvar has a more complex effect on how damages compare with the no-variability damage function. Due to the variability in warming, and for some sectors that have strong seasonal dependence, the present values are smaller than those obtained using the no-variability damage function (Agriculture, Water, Air pollution, Forestry). This highlights the importance not only of incorporating temperature change variability, but also of developing intra-annual, sectoral damage functions. The social cost of carbon (SCC) is a commonly used metric to inform decision-making about the impacts of adding an extra metric ton of CO2 into the atmosphere, or put in a different way, the benefits society can have by mitigating a ton of CO2 24,32,33 . The SCC can be used in cost-benefit analysis of greenhouse gas mitigation projects, and its level can guide setting the optimal (Pigouvian) carbon tax 29 . The SCC was estimated using a simple and physically robust response function of the marginal global mean near-surface temperatures to a one-time pulse of 1 GtC 34,35 released in 2020 (see Methods). This response was added to a set of baseline SSP global temperature scenarios. These scenarios represent a variety of narratives that encompass very high emissions and very high economic growth scenarios (SSP585), one consistent with current trends (SSP370), a middle of the road scenario (SSP245), and two high economic growth scenarios with the successful implementation of stringent mitigation policies (SSP126 and SSP119). Figure 4a shows the estimated SCC values per ton of CO2 for each of the damage functions discussed in this paper and for the SSP585, SSP370, SSP126 and SSP119 scenarios. These estimates are based on a 4% discount rate, but a sensitivity analysis based on 3% and 1.5% discount rate is also included Table S12. In all cases the results show that including spatial variation and temporal variability (S&Tvar) leads to considerably higher SCC values, followed by the S variability damage function, and that temporal variability alone has a marginal effect on SCC. As expected, the largest values of SCC occur under the SSP585 reaching about US$62.24/tCO2 in the case of the no-variability damage function and up to US$78.34/tCO2 with the S&T variability damage function. The S&T estimate for a 3% discount rate is US$120.10/tCO2 (Table S12), which is more than twice that of the US government under the Biden administration for this discount rate 36 , and similar to the price of the carbon tax in Sweden (US$130/tCO2) which is one of the highest in the world. These figures provide further support substantial climate taxes as high as those in countries with the most stringent climate policy 29 . Despite the large changes in climate this trajectory implies, the SCC values for the SSP370 are the lowest due to the much slower economic growth that is associated with such socioeconomic scenario. Using the no-variability damage functions, the SCC value reaches US$25.08/tCO2 compared with US$31.57/tCO2 using the S&T variability damage functions. The SSP2-4.5 and the SSP1-2.6 and SSP1-1.9 imply considerably different levels of warming: while the SSP2-4.5 leads to about 2.72ºC increase in 2100 with respect to preindustrial values, the warming under the SSP1-2.6 and SSP1-1.9 is about 1.77ºC and 1.43ºC, respectively.
At the end of this century, the GDP levels under the SSP2 and SSP1 are very similar, but economic growth is faster during midcentury in the SSP1. Due to the combinations of warming and growth rates, the SCC value is slightly higher under the SSP1-2.6 scenario reaching about US$38.76/tCO2 (US$30.79/tCO2) using the S&Tvar (None) damage functions; for the SSP2-4. The SCC was decomposed into 15 sectors, which add up only approximately to the SCC totals described in the previous paragraph due to differences in warming between the various months and annual warming (Table S13). A sensitivity analysis of the discount rate is presented in Tables S14-S15. Figure 4b shows the results for the SSP5-8.5 scenario obtained using the four sets of monthly sectoral damage functions (results for all SSPs considered are shown in Figure S5). The sector with the highest SCC value is Health (S&Tvar: US$15.96/tCO2 ; None: US$11.89/tCO2) which represents about 20% of the total SCC, followed by Amenity (S&Tvar: US$13.92/tCO2; None: US$10.37/tCO2; about 18% of total SCC), Catastrophe (S&Tvar: US$10.35/tCO2; None: US$8.22/tCO2; about 13% of total SCC) and Water (S&Tvar: US$8.50/tCO2; None: US$7.42/tCO2; about 12% of total SCC), while the only sector with negative SCC value is Time use (S&Tvar: -US$12.06/tCO2; None: -US$8.98/tCO2). These estimates allow to identify that about 60% of the total SCC comes from impacts on four sectors (Health, Amenity, Water, and Ecosystems). Most of the damages occur in the winter and summer, and that damages in 3 of these sectors would exceed 0.1% of global GDP per month during this decade (Figure S4). This suggests prioritizing adaptation in these sectors and seasons is key for reducing climate change impacts. For some sectors, the SCC values obtained from the Svar damage functions are greater than those from the S&Tvar functions (Agriculture, Forestry, Water and Air pollution; Tables S13-S15). In these sectors, the SCC values from the Tvar damage functions are smaller than those obtained from the damage functions that consider no variability (None). Similar to what is discussed in the previous section, these results are produced in sectors that have strong seasonal dependence and in which at least some of the relevant months show lower warming than the annual mean warming. Using modified versions of the damage functions from the RICE model (Methods; Table S16), the global SCC estimates were decomposed into 12 world regions. Figure 4c shows the SCC values for the SSP5-8.5 while Tables S17-S19 show the results for the SSP370, SSP245, SSP126 and SSP119, and three discount rates (4%, 3%, and 1.5%) The damages produced by an extra ton of CO2 are distributed in a highly heterogeneous way, concentrating in the most vulnerable regions. The highest SCC values occur for India (S&Tvar: US$14.28/tCO2; None: US$11.35/tCO2), Africa (S&Tvar: US$14.25/tCO2; None: US$11.32/tCO2), and China (S&Tvar: US$11.55/tCO2; None: US$9.17/tCO2) . These three regions account for about 45% to 51% of the global SCC value, depending on the SSP scenario. In contrast, the SCC values are about one order of magnitude lower for Eurasia (S&Tvar: US$0.99/tCO2; None: US$0.79/tCO2), Russia (S&Tvar: US$1.07/tCO2; None: US$0.85/tCO2) and Japan (S&Tvar: US$1.55/tCO2; None: US$1.23/tCO2). The SCC for the most developed regions (US, WEU, OHI) represents 22% to 25% of the global SCC value. Most IAMs and economic assessments of the impacts of climate change rely on the assumption that changes in the mean of a chosen index (e.g., global mean annual temperature) provide an adequate representation of the effects of this phenomenon on natural and human systems. Although this assumption helps keeping IAMs as simple and tractable as possible, it has two important limitations. First, it fails to take advantage of the advancement that climate modelling has had over the last 30 years in producing spatially explicit projections with better representation of low-and highfrequency variability 37 , and of the development of methods to emulate the output of complex climate models with much simpler and low computing costs 17,38 . Second, as is shown here, this simplification produces nontrivial biases in the estimated costs. Contrary to past studies which have focus on evaluating the effects temporal variability, here we show that omitting spatial variation produces large downward biases on the estimates of the costs of climate change, and that the interaction between temporal variability and spatial variation increases these biases. We show that current IAMs can be easily extended to account for these limitations and provide a more accurate representation of the effects of changes in climate both in the mean and its variability. We develop a set of damage functions that are able to account for temporal and spatial variability in temperature change. Moreover, annual, seasonal and monthly sectoral damage functions are proposed to fully take advantage of seasonal and spatially explicit temperature change projections. Results show that accounting for variability in warming affects the expected level of losses, and not only the risk premium, as proposed in previous publications. Spatial variation and temporal variability together always lead to higher present value of losses and SCC estimates when aggregated across sectors. However, for sectors in which there is a strong seasonal dependence and in which the warming in months that a relevant for the sector is lower than the annual average, economic losses can be smaller than under the assumption of no seasonal variability in warming. We show that traditional damage functions that rely on annual mean global temperature change may severely underestimate the expected losses from climate change. Under the SSP585 scenario, by 2100 global annual losses are projected to be about US$17-US$28 trillion dollars higher when spatial variation and temporal variability are accounted for, in comparison to commonly used damage functions that ignore warming variability. Similarly, the corresponding present value of losses during this century (2020-2100) is about US$47-US$66 trillion dollars higher using the S&Tvar damage functions. It is important to underline that the reported estimates are based on conservative calibrations of global damage functions and a relatively high discount rate is used (4%) for the main estimates. Comparing with previous estimates from Calel et al.
there are differences to point out. Calel e al. 18 report that the expected damages do no change when accounting for this factor, but that the risk premium does by US$32 trillion. Here we show that spatial variation and temporal variability do affect the expected value of damages, increasing losses in the range of $47-$66 trillion or about ½ to ¾ of annual global GDP in 2020. These estimates are much larger than those in Calel et al. 18 although they use one of the damage functions that produce highest damages in the literature 39 . Another important difference is that Calel et al. conclude that accounting for temporal variability has no effect on the social cost of carbon. Our results show that ignoring spatial variation and temporal variability produces an important negative bias on SCC estimates. The SCC estimates are also highly sensitive to the assumptions about warming variability embedded in damage functions. If spatial and temporal variability are included, the SCC values are about 25% higher than under the assumptions in traditional damage functions. Moreover, developing monthly sectoral damage functions allow to understand better how sectors with strong seasonal components can be affected by climate change and, in some cases, to avoid over estimating losses. The damage functions developed in this paper are expected to be particularly useful when applied to climate scenarios in which warming spatial and temporal patterns are highly modified. This type of scenarios include climate catastrophes such as the shutdown of the thermohaline circulation, the Amazon dieback, and the disintegration of Wests Antarctica icesheet. Overly optimistic estimates of the costs and consequences of such events may be largely due to current IAMs ignoring large-scale and profound changes in climatic patterns would bring and concentrating only in the changes these events produce over global mean warming. We use the global monthly gridded dataset of near surface temperature projection of the ensemble mean of the Coupled Model Intercomparison Project Phase 6 (CMIP6; one realization for each of the 37 climate models included). The sample period is 1850-2100. The SSP585 scenario was chosen to maximize the climate signal, although sensitivity analyses are presented for the SSP119, SSP126, SSP245 and SSP370 in the Supplementary Information. We used realizations from different climate models with various ensemble sizes: ACCESS-CM2 (ensemble mean of 3 realizations), CanESM (ensemble mean of 25 realizations), CESM2 (ensemble mean of 25 realizations), INM-CM5-0 (one realization), MPI-ESM1-2-HR (one realization), and ACCESS-ESM1-5 (one realization). All simulations are expressed as anomalies with respect to the 1850-1880 period. For each model, we used two different procedures to separate the climate signal from the weather (irregular) component. The first procedure (SP) is particularly convenient for IAMs. It involves calculating the scaling patterns from general circulation climate models' simulations, which allow to produce spatially explicit climate projections using a simple zero-dimensional climate model 14,16,40 . To generate these scaling patterns, we use ordinary least squares to estimate the slope parameters in regressions between local and global temperature for each grid cell with geographical coordinates 𝑖, 𝑗 41 . To construct spatially explicit future scenarios, these scaling patterns are then scaled with global temperature change projections. All global temperature projections used here were generated with MAGICC7 22,42   . The use of different models to obtain the corresponding spatial patterns of change, scales by a single realization of global temperature (per SSP investigated), allows to evaluate the effects of variability in a separate manner from the uncertainty in climate sensitivity. In the second procedure (RM), we approximate the climate signal using a rolling mean smoothing procedure with a 31-year window for each grid cell. This averaging minimizes the effects of high-to low-frequency oscillations and allows to approximate the underlying climate of the chosen region. This procedure also helps to alleviate some of the problems caused by nonstationarity in climate variables 23,43 . The effect of a pulse of 1GtC on global temperatures is approximated using the stylized response function proposed by Ricke and Caldeira 34 , which provides an option that is physically sound and easy to implement in IAMs for estimating SCC values. It is composed of three-exponential function

, where the units of ∆𝑇 𝑡 are mK/GtC. The SCC estimates presented here are US dollars per tCO2, which is converted from tC by multiplying by 12/44. To represent economic growth for the remaining of this century (2020-2100), we use the estimates of global and for the RICE regions GDP produced by the OECD Env-Growth for the SSP5 narrative, included in the SSP Database 44 . These projections were linearly interpolated from 5-year to annual frequencies. The damage function originally proposed by William Nordhaus to transform global mean annual surface temperature increase from pre-industrial times (∆𝑇 𝑡 ̅̅̅̅̅ ) into economic damages (𝐼 𝑡 ) has been updated over the years, but its main features have remained the same [24][25][26][27] . It consists of a power (=2) function of ∆𝑇 𝑡 ̅̅̅̅̅ , scaled by some parameter values that represent the loss in global welfare for a given temperature change and a calibration temperature value (𝐶). The various versions of this damage function can be succinctly described by the following equation:

Where 𝑎 represents the losses estimated for an increase of 𝐶 in global temperature and 𝛼 = ( 𝑎 𝐶 ). By means of a review of published estimates based on Tol (2009,2014; see Supplementary Information S1), the selected values for 𝑎 and 𝐶 we use here are 2.01% GDP and 2.5ºC, respectively. It is important to note that ∆𝑇 𝑡 ̅̅̅̅̅ denotes the average (expected value) of temperature change over the entire planet and a sufficiently long period of time (e.g., 30 years). Making use of the derivation of the variance, 𝑣𝑎𝑟(𝑋) = 𝐸[𝑋 2 ] -𝐸[𝑋] 2 , so that 𝐸[𝑋 2 ] = 𝐸[𝑋] 2 + 𝑣𝑎𝑟(𝑋), the damage function in equation ( 1) can be written as:

Note that if 𝑣𝑎𝑟(∆𝑇 𝑡 ) = 0, equation ( 2) is reduced to equation (1).
That is, Nordhaus' original formulation implicitly assumes that there is no temporal nor spatial variation in the temperature change. Alternatively, it assumes that the variance of the change in temperature changes in direct proportion to the change in temperature squared: [𝑣𝑎𝑟(∆𝑇 𝑡 )] = 𝛽(∆𝑇 𝑡 ̅̅̅̅̅ ) 2 . To incorporate the effects of variability in temperature, the assumption of zero variance/proportionality is relaxed. The damages that include temperature variability are denoted 𝐼 𝑡 * and are computationally more efficient to calculate as follows:

The difference between equations ( 1) and ( 3) is that the first calculates the damages associated with the average temperature change over space and time domains, while the second represents the average damages of temperature change over both dimensions. As a simple illustration consider the case of two numbers that are not the same, for example 1 and 3. If the average of these two numbers is taken first and then raised to the power of 2, the result is

In contrast, if the numbers are raised to the power of two and then the average of the two is taken, the result is (

If there is variability in the data, the average of a sum of squares will be larger than the average of the same set of numbers risen to the power of 2. Since damage functions evaluate the losses caused by changes in climate, the value of 𝐼 𝑡 * in the reference year (i.e., preindustrial times) is subtracted from the 𝐼 𝑡 * sequence. The damage functions expressed by equations 1-3 assume that the impacts of climate change are uniform over the year. This is expected to be true for some of them such as sea level rise, which is driven by thermal expansion and ice melt, that is cumulative warming rather than warming in a specific month. However, such impacts as the demand for space heating or cooling, or cold-and heat-related mortality and morbidity are tied to the projected warming in winter and summer, respectively. We can thus improve on the naive calibration of Equation (3). We first compiled results from past studies in which the economic impacts of climate change have been estimated (Supplementary Information S1), particularly those in which sectoral decomposition is included. There are four studies that provide estimates for sectoral impacts 19,29,45 , and they are summarized in Table S1. The total damage from all sectors is estimated to be 2.01% of global GDP for a 2.5ºC warming, where 40% and 46% of them are market and non-market damages, respectively, and 13% correspond to risk premium (catastrophes) 29 . As a first step to derive seasonal and monthly damage functions, for each sector seasonal weights are assigned to each sector depending on how active the sector is at that time of year (Table S2). The seasonal damage parameters become:

Where 𝛼 𝑘,𝑗 𝑠 and 𝜔 𝑘,𝑗 𝑠 are the parameter of the seasonal damage function for season and the weight for season 𝑘 and sector 𝑗, and 𝛼 𝑗 is the original damage parameter obtained form the literature for sector 𝑗. Using these weights, the total damages per season for 2.5ºC warming are estimated to be 0.29%, 0.92%, 0.29% and 0.52% for spring, summer, fall and winter, respectively. The seasonal weights and the differences in total GDP per hemisphere are used to derive monthly damage functions as described in the following equation:

Where 𝛼 𝑙,𝑗 𝑚 is the monthly damage coefficient for month 𝑙 and 𝑙 ∈ 𝑘, 𝜔 ℎ 𝐺𝐷𝑃 is the ratio of GDP produced in hemisphere ℎ ($77,800 and $9,470 billion USD for the northern and southern hemispheres, respectively) and global GDP ($87,270 billion USD). The monthly damage coefficients imply economic losses (%GDP) for a 2.5ºC in warming of about 1.10% for December, January and February, 0.55% for March, April and May and, 1.69% for June, July and August, and 0.55% for September, October and November. Note that the sum of monthly (as well as seasonal) damages equals the annual total 2.01%. Figure S1 shows the seasonal and monthly aggregated damages as percent of GDP for 1ºC to 6ºC increase in annual mean global temperatures. We adapt the regional damage functions included in the RICE model 19,20 to include temporal and spatial variability. These 12 regional functions project economic damages based on changes in annual mean global temperature. The regions included in the RICE model are the US, Western Europe (WEU), Japan, Russia, Eurasia, China, India, Middle East (MEAST), Africa, Latin America (LAM), other Asian countries (OASIA), and other high income countries (OHI) 19 . As first step, we standardized the specification of these regional functions, as half of them (China, India, EAST, Africa, LAM, and OASIA) include a linear term in addition to the quadratic function in equation ( 1). We use an optimization method to find the values of parameter 𝛼 𝑟 for each region r that minimizes the squared differences between the original functions and the one that only includes a quadratic term for a range of global temperature change from 0 to 6ºC. This approximation provides very good fit, with 𝑅 2 values larger than 0.99 in all cases. The original and fitted parameter values are provided in Table S16. For each region, the proposed damage functions that include temporal and spatial variability are obtained following equations (1) to (3). Note that results from global damage functions such as the one in (1) and the aggregated total obtained from regional damage functions (like as those in RICE) do not exactly match 19 . To avoid discrepancies between results here we scale regional results to exactly match those based in equation (1). The scaling factor is the ratio of the results from the global damage function and the sum of those obtained for each region with RICE. This ensures that regional and global results are exactly consistent, while it maintains the regional differences generated by the region-specific damage functions. Figures are in billion US$2005.
S&Tvar, Tvar and Svar denote the use of damage functions that include both spatial variation and temporal variability, and temporal variability and spatial variation alone, respectively. The column labeled None shows the results produced with the damage function that does not include any type of variability in temperature change. Results are based on the SP approach. Numbers in parenthesis denote benefits. Calculations are based on a 4% discount rate. Agriculture Figures are in billion US$2005. S&Tvar, Tvar and Svar denote the use of damage functions that include both spatial variation and temporal variability, and temporal variability and spatial variation alone, respectively. The column labeled None shows the results produced with the damage function that does not include any type of variability in temperature change. Results are based on the SP approach. Numbers in parenthesis denote benefits. Calculations are based on a 4% discount rate. Agriculture   Figure S3. Projected economic damages (% GDP) as a function of changes in global temperature according to the sector level monthly damage functions and the SSP585 scenario. S&Tvar (red), Tvar (dotted red) and Svar (green) denote the use of damage functions that include both spatial and temporal variability, and temporal and spatial variability alone, respectively. The slashed blue line labeled None shows the results produced with the damage function that does not include any type of variability in temperature change.