# Climate migration amplifies demographic change and population aging

## Significance

## Abstract

**P**

_{t + 1}is the population at time

*t*+ 1,

**S**

_{t}is a Leslie matrix containing the age–sex-specific probabilities of fertility and survival, and

**M**

_{ty}is a matrix containing the proportion migrating from county to county under

*y*amount of SLR. See

*Methods*and

*SI Appendix*for details.

**P**

_{t}comes from the National Vital Statistics System (NVSS) US Census Populations with Bridged Race Categories Dataset (https://seer.cancer.gov/popdata/download.html).

**S**

_{t}is populated with cohort-change ratios (CCRs) from the NVSS population data. Finally, migration research suggests that “climate migration” is likely to follow preexisting migration flows (28, 29), leveraging embedded social capital and kin networks in destination decision-making (5). Thus,

**M**

_{ty}contains the proportion migrating from each county to each county based on the IRS county-to-county migration data (30) and adjusted to account for the

*h*proportion of each county inundated by

*y*SLR.

**M**

_{ty}reflects the percentage of population we anticipate will be displaced by SLR from a parsimonious displacement model. We arrive at this reduction with a simple, parsimonious model fit with US counties between 1980 and 2019 with verified large (>4

*σ*) population declines (

*n*= 48 such county-years). We use the Spatial Hazard Events and Losses Database for the United States (SHELDUS)—a county-level hazard dataset containing information about the direct losses (property and crop losses, injuries, and fatalities) caused by a hazard event (thunderstorms, hurricanes, floods, wildfires, tornadoes, flash floods, earthquakes, etc.)—to verify these large population declines. We estimate exposure to SLR as a percentage of the population in each county using a bathtub model of inundation based on the population-weighted area under Representative Concentration Pathways (RCPs) 2.6, 4.5, and 8.5 (for 2100 SLR amounts of 0.5 m [0.29–0.82], 0.59 [0.36–0.93], and 0.79 [0.52–1.2] 90th percentile prediction intervals). For simplicity and clarity, we report most results using RCP4.5.

**M**

_{ty}and

**S**

_{t}are projected using ARIMA(0,1,1) (

**S**

_{t}) and ETS (error, trend, seasonal models;

**M**

_{ty}) to capture potential changes in demographic rates.

## Results

^{*}people (Displacement model) to migrate between 2020 and 2100 under all feasible Representative Concentration Pathway-Shared Socioeconomic Pathway (RCP–SSP) scenarios (Figs. 2

*A*and 3). The potential migration of 0.4–10 M people (Amplification model) leads to a demographic amplification of 5.7–53 M people. To illustrate the demographic amplification of SLR, RCP4.5–SSP2 manifests as 1.5 M people migrating to other places for a net demographic change of 1.5 M people (Displacement model). However, accounting for population dynamics under RCP4.5–SSP2 (Amplification model), the actual demographic impact of 1.5 M migrants is 15 M. This is due to two interacting effects: 1) a compounding or “domino effect” where destination counties attract more people over time and origin counties attract less and the interaction of migrants with the other demographic component processes of fertility and mortality. 2) This total demographic effect or amplification is considerably more pronounced than the simple displacement effect—10 times larger under RCP4.5–SSP2 and between 5.3 and 18 times larger under all feasible RCP–SSP combinations. Under all plausible RCP–SSP scenarios except RCP4.5–SSP3, SLR leads to a demographic amplification of at least 10 million persons (Table 1).

RCP | SSP1 | SSP2 | SSP3 | SSP4 | SSP5 | |
---|---|---|---|---|---|---|

8.5 | SLR amount (meters) | 0.79 [0.52–1.2] | ||||

Migrants (millions) | 3.4 [1.3–10] | |||||

Demographic amplification (millions) | 28 [17–53] | |||||

4.5 | SLR amount (meters) | 0.59 [0.36–0.93] | ||||

Migrants (millions) | 1.5 [0.65–4.2] | 1.5 [0.63–4.1] | 0.84 [0.36–2.3] | 1.2 [0.51–3.3] | 2.2 [0.96–6.2] | |

Demographic amplification (millions) | 15 [10–27] | 15 [10–26] | 8.6 [5.7–15] | 12 [8–21] | 23 [15–41] | |

2.6 | SLR amount (meters) | 0.5 [0.29–0.82] | ||||

Migrants (millions) | 1.2 [0.55–3.5] | 1.2 [0.54–3.4] | 0.96 [0.44–2.8] | 1.8 [0.82–5.1] | ||

Demographic amplification (millions) | 14 [9.7–24] | 14 [9.5–24] | 11 [7.6–19] | 21 [15–37] |

Here, we compare the expected number of migrants with the demographic amplification for all reasonable RCP–SSP combinations. All numbers in parentheses are the 90th percentile prediction interval. Regardless of the RCP-SSP combination, matrix population models suggest that millions of people will shift in the United States due to sea-level rise.

*E*and

*F*are counties experiencing significant population declines over the entire time horizon (vulnerable counties). Our results for Miami-Dade, FL, under RCP4.5 suggest 28.3 K [4.4–204.2 K] migrants by 2100 but a total demographic impact of 243.9 K [72.9 K–1.1 M] fewer people. We see similar results for Dare, NC, as well (8.5 K [1.5–22.9 K] migrants and a 39.8 K [14.6–70.0 K] total demographic impact). Relatively small displacements can cascade into considerably larger demographic changes.

*G*and

*H*are counties close to heavily threatened areas. We see two separate population trajectories and demonstrate the amplification of population change. McIntosh, GA, is initially a recipient county, as people from nearby heavily affected areas migrate into nearby counties. But as the century wears on, this recipient county turns into a vulnerable county, and out-migration begins to turn into negative population growth. The difference between the simple, Displacement model and the demographically integrated Amplification model suggests markedly different population growth trajectories. The Displacement model suggests that SLR migration in McIntosh will yield more population than the climate agnostic Base projection into the end of the century. Conversely, the Amplification model has the population greatly increasing before population growth turns negative after 2070.

*B*–

*D*are examples of climate destinations (counties considered “climate havens”). We see a similar, though reversed, pattern of population trajectories to vulnerable counties. Here, Rutherford, TN (245.1 K [78.5–852.3 K] amplified demographic change versus 34.9 K [8.2–197.5 K] displaced migrants), Douglas, CO (377.7 K [117.8–1.6 M] vs. 32.3 K [5.2–230.3 K]), and Washington, OR (131.9 K [50.9–377.8 K] vs. 12.9 K [3.8–55.1 K]), exhibit considerably larger demographic impacts when accounting for population dynamics compared to just the Displacement model. These results echo suggestions of emerging climate havens in “safer” areas (19). If cities prepare solely for climate migrants, they are likely to underprepare for the demographic change associated with climate in-migration.

## Discussion

**2**) requires minimal translation of climate hazards into age-specific displacement. Future scientists could implement extreme heat/humidity combinations, tropical cyclones, or water availability as either singular climate hazards or in a multihazard model. To our knowledge, few, if any, demographic projection models account for climate change impacts. This type of modeling approach has great potential to better understand the integration of population processes and climate change.

## Materials and Methods

**1**is the general form for our projections.

### Parsimonious, One-Dimensional, Age-Specific Displacement Model.

*t*-static are considered “true” outliers (|

*τ*|≥4;

*P*-value < 0.000063). We chose this threshold to minimize the probability of committing a type I error (or claiming that an outlier is true when it is, in fact, not).

*σ*or greater. We then use the SHELDUS (48), a county-level hazard dataset for the United States which contains information about the direct losses (property and crop losses, injuries, and fatalities) caused by a hazard event (thunderstorms, hurricanes, floods, wildfires, tornadoes, flash floods, earthquakes, etc.) for the period 1960 to the present. We use SHELDUS to ensure that the county time periods we identify as statistical outliers with population losses experienced an environmental hazard in that county-year with per capita hazardous losses in excess of the 50th percentile. This is to ensure that the outlier population losses that we detect are associated with a hazard rather than other forces, such as economic forces.

_{n}

*P*

_{x, t}is the population aged

*x*to

*x*+

*n*in time

*t*, and

_{n}

*P*

_{x − k, t}is the population aged

*x*−

*k*to

*x*+

*n*−

*k*in time

*t*, where

*k*refers to the time difference between time periods. Since mortality must decrement a population, any CCR above 1.0 implies a net-migration rate in excess of the mortality rate and a growing population. There is special consideration for both the initial age group without a preceding age group and the final, open-ended age interval without a proceeding age group (see (49) for additional details).

*x*,

*Δ*

*C*

*C*

*R*

_{x}=

*C*

*C*

*R*

_{x, t}/

*C*

*C*

*R*

_{x, t − 1}, and the percentage decline in the total population compared to the counterfactual in the outlier detection method, $\mathrm{\Delta}{P}_{t}=\widehat{{P}_{t}}/{P}_{t}$:

*h*is the log(

*Δ*

*P*

_{t}) and shows a quadratic relationship with the logarithm of the change in CCRs by age.

*x*refers to five-year age groups: 0–4, 5–9,…, 85+. This is a similar model and approach to Wilmoth et al.’s (51) flexible, one-dimensional mortality model based on the similarity between age-specific mortality rates and infant mortality.

*SI Appendix*, Table S1 depicts the r-squared values between log(

*Δ*

*C*

*C*

*R*

_{x}) and log(

*Δ*

*P*

_{t}). The age groups with the lowest r-squared values are young adult males aged 20–39 and those in the open-ended age interval (80+), suggesting that these age/sex groups react to environmental signals the least predictably.

*h*

_{ty}=

*l*

*o*

*g*(

*Δ*

*P*

_{ty}), we can estimate age-specific changes in CCRs after an environmental event by simply applying the following formula:

*C*

*C*

*R*

_{xt}based on the log(

*Δ*

*P*

_{ty}) under SLR amount

*y*. In this case, we drop the intercept (

*a*

_{x}) from the estimation procedure to ensure that a 0% decline in population yields a corresponding 0% change in the CCR. Multiplying the result from the model with the CCR in the year prior yields the anticipated change in the CCR. These changes in CCRs can then be applied to any time series of population values to generate an anticipated population.

*SI Appendix*, Figs. S1 and S2, and

*SI Appendix*, Table S1 show the accuracy of our fitted, one-dimensional model. Regarding the total population in our 48 counties, our model performs well with an

*r*

^{2}value of 0.996 and performs well regardless of population size. Regarding each individual

*P*

_{x}group, our model still performs quite well with an

*r*

^{2}of 0.995. Just like with the total population, the accuracy of our model does not depend on population size.

### Sea-Level Rise Exposure.

*h*in Eq.

**2**, we employ inundation modeling (25), which assumes that people who are underwater 100% of the time must migrate. We estimate these populations using airborne lidar-derived digital elevation models (DEMs) produced by National Oceanic and Atmospheric Administration (NOAA) and supplemented with both the US Geological Survey (USGS) Northern Gulf of Mexico Topobathymetric DEM in Louisiana and the USGS National Elevation Dataset in the fraction of land not covered by other sources (see ref. 52 for details on the construction of the DEMs).

*n*= 81,815) located in coastal counties (

*n*= 406) expected to experience any probability of flooding. We use probabilistic SLR projections (35, 38) that are closely aligned with the IPCC for our water heights. To calculate the land area under a given water height, we simply threshold the DEM to find pixels below

*S*

*L*

*R*

_{yt}, where

*y*is the projected height of SLR in a given year

*t*. For each CBG, we simply calculate the percentage of its pixels on dry land (defined in the National Wetland Inventory (53)) covered by the inundation surface and multiply this percentage by the total population in the CBG in year

*t*to produce the total number of people at risk to a given amount of SLR in a given year. We then aggregate these CBGs to the county level to calculate the percentage of people in a given county at risk of inundation.

*Δ*

*P*

_{t}from our Eq.

**2**above, where

*Δ*

*P*

_{t}=

*P*

_{t, y}/

*P*

_{t}. In this case,

*P*

_{t, y}/

*P*

_{t}is the percentage of the population in any county at risk of inundation under a given water height

*y*. Such a calculation allows us to seamlessly combine our parsimonious Displacement model with our matrix population model.

### Matrix Population Models.

**P**

_{t}refers to the population matrix containing

*x*age groups and

**S**

_{t}contains the age-specific fertility (

*F*) and mortality rates (

*S*), where

*S*

_{x}refers to the survival probabilities for age group

*x*and

*F*

_{x}refers to the fertility rates.

**P**

_{t}is a vector of length 113,148, and

**S**

_{t}and

**M**

_{ty}are 113,148 × 113,148 matrices. We use CCRs to populate our

*S*

_{x}values in each matrix and child–woman ratios to populate our

*F*

_{x}values. The values come from the NVSS Bridged Race Categories dataset. To project the CCRs, we employ an ARIMA model for forecasting equally spaced univariate time series data. We use an ARIMA(0,1,1) model, which produces forecasts equivalent to simple exponential smoothing. All projections were undertaken in R (44) using the forecast package (46).

*c*, age group

*x*, sex

*s*in an ARIMA(0,1,1) to create

*C*

*C*

*R*

_{cxst}for time periods

*t*+ 1. The initial

**P**

_{1}and

**S**

_{t}matrices use the 2019 NVSS data.

**M**

_{ty}matrix. These data come from the Internal Revenue Service (IRS) county-to-county migration files for 1990–2018 (see refs. (54) and (55) for descriptions of this data). We populate

**M**

_{ty}as ${M}_{t}\xb7(1-{e}^{{\widehat{\beta}}_{x}{h}_{\mathit{ty}}+{\widehat{c}}_{x}{h}_{\mathit{ty}}^{2}})$ and where

*i*=

*j*as ${e}^{{\widehat{\beta}}_{x}{h}_{\mathit{ty}}+{\widehat{c}}_{x}{h}_{\mathit{ty}}^{2}}$.

*M*

_{t}is a vector containing the proportion migrating from county

*i*to county

*j*in the IRS migration data, thus ensuring the columns of

**M**

_{ty}sum to 1.0.

^{17}= 5x, 1.3

^{17}= 86x). Projecting the migration system itself is not subject to such exponential drift.

**M**matrix above where nonmigrants (i.e., those migrating from

*i*→

*i*) are included. The result is the fraction of individuals surviving from age group

*x*who migrate from

*i*→

*j*.

*SI Appendix*.

## Data, Materials, and Software Availability

## Acknowledgments

### Author contributions

### Competing interests

## Supporting Information

- Download
- 330.95 KB

## References

*Climate change: the 1990 and 1992 IPCC assessments, IPCC first assessment report overview and policymaker summaries and 1992 IPPC supplement*(IPCC, Geneve, 1992).

*Groundswell*(The World Bank Group, London, 2018).

*Proc. Natl. Acad. Sci. U.S.A.*

**117**, 11350–11355 (2020).

*Report on the Impact of Climate Change on Migration*(US White House, Washington DC, 2021).

*Nature*

**478**, 447–449 (2011).

*Nat. Clim. Change*

**7**, 321–325 (2017).

*Environ. Res. Lett.*

**13**, 1 (2018).

*Earth’s Future*

**9**, 1 (2021).

*Popul. Space Place*

**21**, 54–67 (2015).

*Demography*

**25**, 355–370 (1988).

*Glob. Environ. Change*

**21**, S94–S107 (2011).

*Glob. Environ. Change*

**21**, S70–S81 (2011).

*Int. Migr.*

**49**, e224–e242 (2011).

*Sustain. Sci.*

**9**, 331–345 (2014).

*Monthly Lab. Rev.*

**131**, 32 (2008).

*Popul. Res. Policy Rev.*

**1**, 1 (2021).

*Popul. Environ.*

**33**, 28–54 (2011).

*Nat. Rev. Earth Environ.*

**1**, 28–39 (2020).

*J. Environ. Stud. Sci.*

**11**, 465–480 (2021).

*Glob. Environ. Change*

**42**, 181–192 (2017).

*et al*., IPUMS USA: Version 11.0 (2021) Version Number: 11.0 Type: dataset.

*Nature*

**271**, 321 (1978).

*Environ. Urban.*

**19**, 17–37 (2007).

*Oceanography*

**24**, 144–157 (2011).

*Proc. Natl. Acad. Sci. U.S.A.*

**112**, 13508–13513 (2015).

*PLOS ONE*

**10**, e0118571. (2015).

*PloS One*

**15**, e0227436. (2020).

*Glob. Environ. Change*

**21**, S50–S58 (2011).

*Demography*

**57**, 1437–1457 (2020).

*J. Econ. Perspect.*

**25**, 173–196 (2011).

*Geosci. Model Dev.*

**9**, 3461–3482 (2016).

*Clim. Policy Initiative*

**32**, 1–38 (2014).

*J. Environ. Plann. Manag.*

**63**, 2102–2143 (2020).

*Nat. Clim. Change*

**6**, 691–695 (2016).

*Earth’s Future*

**2**, 383–406 (2014).

*Science*

**315**, 368–370 (2007).

*Proc. Natl. Acad. Sci. U.S.A.*

**116**, 11195–11200 (2019).

*et al*., Global and regional sea level rise scenarios for the United States (2017).

*Science*

**372**, 1279–1283 (2021).

*Wiley Interdiscip. Rev.: Clim. Change*

**10**, e565 (2019).

*Earth’s Future*

**2**, 383–406 (2014).

*Demography*

**3**, 537–544 (1966).

*J. Am. Stat. Assoc.*

**88**, 284–297 (1993).

*R: A Language and Environment for Statistical Computing*(R Foundation for Statistical Computing, Vienna, Austria, 2019).

*Tsoutliers: Detection of Outliers in Time Series*. R package version 0.6-8 (2019).

*et al*.,

*Forecast: Forecasting functions for time series and linear models*. R package version 8.10 (2019).

*J. Stat. Softw.*

**26**, 1–22 (2008).

*The Spatial Hazard Events and Losses Database for the United States*. Version 17.0. [Online Database]. Phoenix, AZ: Center for Emergency Management and Homeland Security, Arizona State University (2018).

*Sci. Data*

**6**, 1–15 (2019).

*Popul. Res. Policy Rev.*

**29**, 47–63 (2010).

*Popul. Stud.*

**66**, 1–28 (2012).

*Nat. Commun.*

**10**, 1–12 (2019).

*National wetlands inventory website. U.S. department of the interior, fish and wildlife service*, Washington, D.C. http://www.fws.gov/wetlands/ (2012).

*Demogr. Res.*

**40**, 1153–1166 (2019).

*Popul. Res. Policy Rev.*

**41**, 1–12 (2021).

*Forecasting with Exponential Smoothing: The State Space Approach*(Springer Science & Business Media, 2008).

*mathewhauer/SLR-mig-proj: Verion 1.1*(2022).

## Information & Authors

### Information

#### Published in

#### Classifications

#### Copyright

#### Data, Materials, and Software Availability

#### Submission history

**Received**: April 11, 2022

**Accepted**: October 16, 2022

**Published online**: January 8, 2024

**Published in issue**: January 16, 2024

#### Keywords

#### Acknowledgments

##### Author Contributions

##### Competing Interests

#### Notes

### Authors

## Metrics & Citations

### Metrics

#### Citation statements

#### Altmetrics

### Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

## View Options

### View options

#### PDF format

Download this article as a PDF file

DOWNLOAD PDF### Get Access

#### Login options

Check if you have access through your login credentials or your institution to get full access on this article.

Personal login Institutional Login#### Recommend to a librarian

Recommend PNAS to a Librarian#### Purchase options

Purchase this article to access the full text.