Modeling Human Migration Under Environmental Change: A Case Study of the Effect of Sea Level Rise in Bangladesh

Sea level rise (SLR) could have catastrophic consequences worldwide. More than 600 million people currently living in coastal areas may see their livelihood at risk and choose to migrate in the near future. Predicting when, how, and where people could migrate under environmental change is critical to devise effective policy initiatives and improve our preparedness. Here, we propose a modeling framework to predict the effect of SLR on migration patterns from easily accessible geographic and demographic data. The framework adapts the radiation model to capture unwillingness or inability to migrate of affected residents, as well as return migration and cascading effects in migration patterns. We apply the mathematical model to study internal migration in Bangladesh, where we predict a complex and counterintuitive landscape of migration patterns between districts. Our predictions indicate that the impact of SLR on 816,000 people by 2050 will trigger cascading effects in migration patterns throughout the entire country. The population of each of the 64 districts will change, leading to a total variation of 1.3 million people. Migration from inundated regions in the center will trigger non‐trivial patterns, including a reduction in the population of the district of the capital Dhaka.

Although there is large variability in future sea level projections, due, for instance, to the uncertainty in anthropogenic emissions, there is consensus on the potentially catastrophic worldwide impact of SLR (Kopp et al., 2019). Coastal regions in the United States might experience migrations "of a magnitude similar to the 20th century Great Migration of southern African-Americans" (Hauer et al., 2016). Similarly, developing countries, like Bangladesh, may see the majority of their population experience SLR-related impact by the middle of the century, potentially triggering massive migrations (Neumann et al., 2015). Not only is formulating reliable predictions of environmental migration critical for devising effective policy initiatives around the implications of SLR, but also it plays a fundamental role in improving our preparedness for future migration patterns.
Several mathematical models have been proposed to predict the impact of SLR on migration patterns (Hauer et al., 2019). For example, Hassani-Mahmooei and Parris (2012) established an agent-based model to examine environmental migration in Bangladesh. Their results suggest a migration in 2050 toward east and northeast districts, along with a continuous absorption of migrants in current urban areas. Through a linear probability model, Chen and Mueller (2018), determined that gradual increases in soil salinity due to coastal flooding could be a stronger factor for environmental migration in Bangladesh than inundation alone. Agent-based modeling was also utilized by Karanci et al. (2017) to predict housing decisions in the coastal town of Nags Head in North Carolina, USA. While having considerable technical merit, these modeling approaches are relatively complex to implement, with several parameters to be empirically calibrated toward accurate predictions.
The study by Davis et al. (2018) makes an important contribution to the literature, by proposing a universal model of human migration under environmental change. The model constitutes a valuable trade-off between complexity and accuracy to predict the influence of SLR on migration fluxes. The approach adapts the radiation model by Simini et al. (2012), which has been shown to accurately describe mobility patterns as a function of population distribution. Specifically, the radiation model predicts that the flux between two regions is controlled by their populations and those of all the neighboring regions. The adaptation by Davis et al. (2018) reduces the value of the populations in the expression of the fluxes to account for human migration triggered by SLR. The model is parameter-free, whereby it only takes as input the population distribution and the population in inundated regions for the prediction of migration fluxes. The approach was illustrated through the study of future environmental migration driven by SLR in Bangladesh.
Albeit offering the first parameter-free approach for elucidating environmental migration, the original approach presented by Davis et al. (2018) has a number of limitations that are, in fact, discussed by the authors in their manuscript. In its present incarnation, the model does not allow for considering unwillingness and inability to migrate, which is widely documented in the technical literature (Hauer et al., 2019). For example, previous studies on Bangladesh have indicated that residents of affected regions display a complex coping response to hazard, which may not trigger migration (Barman et al., 2012;Hutton & Haque, 2003). Should one be willing to migrate, the route to migration is still hampered by socioeconomic factors, where the poorest and least mobile will be "trapped populations" (Black, Bennett, et al., 2011).
In addition, the model neither accounts for cascading effects in environmental migration nor for return migration (Hauer et al., 2019). As large populations are displaced due to environmental change to a new region, they may conflict with existing residents who may see their economic resource base threatened (Reuveny, 2007). Such a conflict may, in turn, trigger consequent migrations in a cascading effect. Displaced people may also express a desire to return to their homes, as seen for example, among New Orleans' residents who returned to New Orleans after Hurricane Katrina for a sense of place (Chamlee-Wright & Storr, 2009). Both cascading effects and return migration will ultimately contribute to a richer landscape of migration patterns, in which local effects of environmental change in any region will have global consequences.
Here, we propose a mathematical model to study migration patterns due to environmental change that seeks to overcome these limitations. Each region is treated as a node in a weighted network (Estrada, 2012), whose directed links represent migrations between regions. We encapsulate demographic and geographic data in each node through the population and position of the corresponding region. The population of a region may increase or decrease due to incoming or outgoing links, respectively. The migration rate along each link is controlled by the population of the endpoints, as well as the population of any other regions that is DE LELLIS ET AL. 10.1029/2020EF001931 2 of 14 geographically closer to the origin than the destination, following the radiation model (Simini et al., 2012). The migration flux is then computed by multiplying the corresponding rate by the population of the origin that is willing and able to migrate under the effect of environmental change.
Different from the approach of Davis et al. (2018), the model contains one parameter that captures baseline migration between regions, in terms of the fraction of lifetime migrants who are born in a region and live their lives in another region. The presence of environmental change is modeled through a localized shock onto this baseline parameter. The more severe is the shock, the larger is the fraction of people that could be displaced in a region, thereby changing the global landscape of migration patterns. Adapting techniques from the study of vulnerability and leadership in network dynamical systems (Fitch & Leonard, 2015;Pagnier & Jacquod, 2019;Porfiri & Frasca, 2018), we propose a sensitivity analysis to examine the specific effect of a local shock on global migration patterns. In this context, we apply a perturbation at a given node (that is, an increase in the number of individuals who are willing and able to migrate) and study its cascading effects in any other node. We introduce a resilience index to measure the extent to which a localized perturbation induces migration everywhere in the network.
By linking the results of the resilience analysis with predictions about the spatial distribution of environmental change, we quantify the extent of the potential threat on migration patterns and offer long-term predictions. For example, the same average SLR for a given country may trigger widely different migration patterns, depending on how the population is distributed between coastal and inland regions and between rural and densely populated urban areas. We focus the analysis on the effect of SLR on internal migrations in Bangladesh, demonstrating complex inter-dependencies between SLR distribution and migration patterns between different regions.

Formulation of the Migration Model
Our modeling approach begins by partitioning the overall geographical area where environmental change and the resulting migration take place into N regions, corresponding to the network nodes. At each node, we measure the population at discrete time instants separated by a unitary time step; the population in the i-th region at time k is denoted as n i (k). By assuming that birth and death rate balance each other and that there is not an influx from or outflux toward other areas, we can write the following continuity equation to predict the time-evolution of the population in the network: where J(i → j, k) is the migration flux from origin i to destination j in the time interval [k, k + 1).
To complete the system of equations and afford predictions of environmental migrations, we must relate the migration fluxes to the populations and some salient parameter of environmental change through constitutive equations. We propose the following constitutive equation of the generic flux: Here,  m ( ) ( ) i i n k n k is the number of people living in region i at time k that could migrate in the time interval [k, k + 1) (superscript "m" stands for "mobile"); n(k) is a vector that aggregates all the populations at time k; and r ij (n(k)) ∈ [0, 1] is the fraction of the mobile population of origin i which chooses j as the destination, which we name the migration rate.
We measure the mobile population in each of the N regions with respect to the initial population at time 0, that is, we introduce N parameters The parameters α 1 , …, α N ∈ (0, 1) are used to capture the effect of environmental change, whereby affected regions will register an increase in the fraction of inhabitants that will try to move, that is, a reduction in the corresponding parameter values. In our modeling approach, environmental change causes an increase in the number of people who are willing and able to move: whether or not the migration will take place depends on demographic and geographic data that involved the entire network. To minimize the number of fitting parameters, we hypothesize that in the absence of environmental change α 1 = ⋯ = α N = : α 0 , where α 0 is the baseline value.
Demographic and geographic data contribute to the value of all the r ij (n(k))-s for i, j = 1, …, N. We take r ij (n(k)) to be equal to the probability of migration obtained from the radiation model (Simini et al., 2012), which relates migration tendency from an origin to a destination to the populations of both these regions as well as the population of any neighboring region. Specifically, we impose r n k n k n k n k n k n k n k n ij i j is the set of all the regions which are closer to the origin than the destination, with d ij being the distance between the centroids of districts i and j. In the original radiation model, Equation 4 is derived under the premises that an individual will select its migration destination on the basis of better life opportunities, while seeking to remain close to the origin. In this vein, the summation in the denominator of Equation 4 over the neighbors accounts for intervening opportunities that should be weighted during the decision process.
The structure of the model equations leads to the following three properties, which are easy to verify: 1. The total population is constant, that is, This property follows directly from the continuity equation in Equation 1 by summing over i.
2. For all i = 1, …, N, and for all times k This property follows from the fact that the outflux at time k can never exceed m ( ) i n k .
3. The probability of not migrating at time k can be computed from the radiation model as We utilize the migration model to predict the populations in a network at the steady-state, starting from given initial conditions for the populations n 1 (0), …, n N (0), a chosen distribution of α 1 , …, α N (encompassing changes with respect to the baseline value α 0 ), and prescribed geographical distances between each pair of nodes d ij with i, j = 1, …, N. This step is simply executed by running the recursion Equation 1 with constitutive Equation 2, specified through Equations 3 and 4. By comparing the steady-state populations with the original populations, we assess the extent of net migration within the network.
We note that r ij (n(k)) depends on demographic and geographic data, without explicitly accounting for the size of the time step, whether it is 1 year or 5 years, for example. The derivation by Simini et al. (2012) does not explicitly include time, such that probability of migration is intended with respect to the time scale of the migration process itself. Should the decision process by an individual unfold on the time scale of 5 years, this number should be intended as the probability of migrating within 5 years and a simulation conducted with a resolution of 1 year should artificially reduce the probability of migration by a factor of five. In our simulations, this aspect has a secondary role, whereby we largely focus on the equilibrium where influxes and outfluxes are completely balanced, that is, k → +∞. The same predictions would be obtained by directly n k n and considering Equation 5, that is, by setting the left-handside to zero.

Analysis of the Migration Model
To facilitate the study of the proposed mathematical model of migration, we introduce the matrix function Q(⋅) with zero-column sum, whose ij-th element is where a is an arbitrary vector, with components a 1 , …, a N . By accounting for Equation 7, we write the model equations as where I N is the identity matrix, diag(⋅) is the diagonal form of a vector, and α collates all the values α 1 , …, α N .
As a first step in the analysis of migration patterns, we perform a sensitivity analysis by increasing the number of mobile individuals in the generic region i by one with respect to the baseline value α 0 , such that  where i e is the vector of all zeros except of a one at position i. The perturbation will modify the equilibrium 0 n corresponding to where ν (i) collates the (zero-sum) changes in the population at the equilibrium for all the nodes in the network due to the perturbation at node i.

By linearizing Equation 9, we establish
where    ( , ) Q is a matrix function with zero-column sum, whose generic element is computed as for any pairs of vectors a and b. The expression of Equation 11 can be explicitly obtained by carrying out the derivatives.
We can compute the steady-state value  ( ) i of Equation 10, either as the limit k → +∞ or simply solving for ν(k + 1) = ν(k) with a zero-sum perturbation. By studying  ( ) i , we establish how a small increase in the mobile population of a region may trigger complex migration patterns throughout the network. As a result, one of the regions will see an increase in the number of incoming migrants and another may register a decrease in the population. We quantify the extent of these effects throughout the network by summing the absolute values of all the variations in the population to obtain  ( ) DE LELLIS ET AL.
10.1029/2020EF001931 5 of 14 By comparing the values of  ( ) 1 i ‖ ‖ for different indices i, we can evaluate resilience to a localized shock anywhere in the network. The higher this number is, the stronger will be the global effect caused by a localized shock on migration patterns. To present results in a normalized manner, we calculate the "resilience" to environmental change at region i as Such a sensitivity analysis offers insight into the relative role of each network node on migration patterns, under the premise that only a small fraction of the population will be affected by a shock. To overcome this approximation, we should simulate the complete nonlinear model (Equation 9) with a large perturbation.
Due to the nonlinearity of the problem, the superposition principle does not apply and we should consider at once all the perturbations that are applied throughout the network. Hence, the nonlinear analysis requires to increase the number of mobile individuals as n m (k) = n(k) − α 0 n(0) + n m−ec , where the vector n m−ec collates the estimated number of people that could become environmental migrants at each node (superscript "m-ec" identifies "mobile" population due to the "environmental change"). Finally, by subtracting 0 n from the resulting steady-state solution n, we can estimate the (zero-sum) variations in the population due to environmental migration.

Demographic and Inundation Data in Bangladesh
The administration of Bangladesh is organized in two main levels: The impact of SLR on Bangladesh is potentially disruptive, since 41% of the total population lives in areas where the elevation is lower than 10 m (Neumann et al., 2015). To quantify the population in each district that could be affected by SLR in 2050, we use the estimations made by Davis et al. (2018). , to finally gather the population that is expected to be inundated in each district. We collate these values in the vector n m-ec with N = 64 entries, of which 21 are different from zero in correspondence to potential environmental migration from inundated districts. From a modeling perspective, this vector consolidates the entire effect of environmental change, thereby serving as the sole input to predict migration patterns.

Calibration of the Migration Model on 2011 Division-Level Migration Data
As a baseline for studying the effect of environmental change, we calibrate the model with respect to division-level migration data from (Bangladesh Bureau of Statistics 2015, 2011 Population and Housing Census, 2011) (N = 7). We hypothesize a common value for α 1 , …, α N equal to α 0 , which we identify from the data. We set the initial population in each of the seven divisions from the 2011 Census data and use geographic distances between the centroids of divisions computed from Davis et al. (2018), to implement the . By varying α 0 from 0 to 1, in steps of 0.001, we determine the value that maximizes the R 2 value between the pairs predictions/data and the bisectrix (R 2 = 0.739, p < 0.001), see Figure 1. While it might be preferable to calibrate at a finer spatial resolution and over a time window, data availability challenged pursuing a different approach. A similar analysis has been, in fact, reported by Davis et al. (2018), although with the different goal of validating the fluxes in the radiation model against lifetime migrants, before examining the effect of environmental change.

Results and Discussion
We articulate the assessment of the severity of the potential threat of SLR on Bangladesh in three sequential steps.
First, we apply the sensitivity analysis to district-level data from 2011 Census (N = 64), thereby computing the resilience in Equation 12. To delve into the mechanisms underlying a differential resilience of the country to environmental change, we perform a correlation analysis linking the resilience to a localized shock in a district with the population in its neighboring districts. Neighbors are defined to be those districts whose centroids are within a distance that is less than the average distance between any two districts (207 km). A positive correlation would indicate that highly populated regions could shield migration from a neighboring district, whereas a negative correlation would reveal cascade effects in migration patterns.
Second, we correlate the values of the resilience indices with the corresponding affected populations n m−ec in the inundated districts (N = 21) by 2050 (estimated by Davis et al., 2018). Through this analysis, we seek to evaluate the overall resilience of Bangladesh to SLR. A positive correlation would indicate the more favorable scenario in which inundated regions have a more peripheral role on migration patterns, such that people who will be placed on the move will not trigger cascading migrations throughout the country. On the other hand, a negative correlation would point at a more dramatic scenario, in which SLR will impact critical districts with respect to global migration patterns.
Third, we perform a full nonlinear analysis of migration patterns that accounts for the simultaneous SLR of the 21 inundated districts. Toward this aim, we estimate the steady-state population in all the 64 districts, population 0 n that are predicted from the same initial conditions, but in the absence of the potential environmental migrants n m−ec .

The Location of Environmental Change Will Have a Differential Effect on Migration Patterns in Bangladesh
We begin by examining the resilience of Bangladesh to localized shocks in each of the 64 districts. The map in Figure 2 suggests that shocks in districts in the north, east, and south will have a lesser impact on migration patterns than central and western districts. However, we should acknowledge a rich spatial dependence, with some of the most critical districts that neighbor less influential regions. The highest value of resilience is registered in the capital district of Dhaka in the center of the country, while the lowest is in the central-western district of Narail. Notably, these rankings are robust with respect to a simplistic inclusion of external migration, in the form of an additional, virtual district that is equally accessible from all the 64 districts (see Supporting Information S1 and Figure S1). Figure 2 indicates that resilience to environmental change in a district is inversely related to the population of neighboring district (R 2 = 0.340, p < 0.001). The more a district is surrounded by highly populated districts, the stronger are the resulting migration patterns within the country. This suggests the possibility to mathematically anticipate cascading effects, in which people affected by the environmental change will first migrate to neighboring districts, thereby triggering further migration to the rest of the country. These cascading effects have been documented in the literature (Hauer et al., 2019) and attributed to potential conflict between migrants and existing residents, who may experience lower quality of life due to the incoming migration.

The correlation analysis in
The correlation analysis in Figure 2 utilizes a common value α 0 , but equivalent predictions would be obtained by accounting for small heterogeneity in the values of α 1 , …, α N with respect to α 0 , as discussed in Supporting Information S2 and Figure S2. In particular, variations of up to 10% in α 1 , …, α N never confound the significance of the correlation and the resulting changes in the slope of the fit are within 5% of the nominal value for a common α 0 .

The Geography and Demographics of Bangladesh Could Magnify the Threat Posed by SLR on Internal Migration
The map in Figure 3 illustrates the number of people who will be impacted by SLR in 2050. Districts in the south will be mostly impacted by SLR, with the largest fraction in the central-southern districts of Shariatpur, Munshiganj, and Narayanganj, where over 335,000 people will be affected.
Although the districts that are expected to be impacted by SLR are not those that yield the lowest resilience from Figure 2, their role is expected to be dramatic. From the correlation analysis in Figure 3, we evidence a negative correlation between resilience to localized shocks and the affected population (R 2 = 0.223, p = 0.031). SLR is predicted to affect large segments of the population living in critical districts, which will trigger global migration patterns throughout the country.
The extent of the effect of SLR can be visualized in Figure 4. Therein, we demonstrate that a unitary increase in the mobile population of Shariatpur (the most affected district by SLR in 2050) will trigger mobility in three northern districts toward western districts. The capital of Dhaka should register an outflux of people DE LELLIS ET AL.   due to the increased mobility of the population in Shariatpur. This possibility might appear counterintuitive, based on predictions of previous models (Chen & Mueller, 2018;Davis et al., 2018;Hassani-Mahmooei & Parris, 2012) and empirical evidence by Kartiki (2011), suggesting that migrants prefer to migrate to urban areas.
The basis of this surprising result is related to the gradual effects of SLR, which drive the system to a new equilibrium. Specifically, our analysis focuses on steady-state migration patterns that will emerge in response to the increased tendency of people to move from inundated areas. While in a single time-step, highly populated districts will tend to attract more migrants, as time evolves this incoming population might tend to diffuse to neighboring regions, until reaching the new equilibrium. For example, a one time-step prediction in Equation 1 about the effect of the shock in Shariatpur (corresponding to a unitary increase in the mobile population) will consist of a normalized increase of 0.048 on Dhaka's population, in contrast with the steady-state decrease of 0.224 shown in Figure 4.
The same phenomenon is not observed in response to environmental change in the peripheral region of Cox's Bazar shown in Figure 4, where more than 23,000 people are estimated to be impacted by SLR. In this case, the neighboring districts are not expected to see their population change, while the central districts of Dhaka, Rajbari, and Narsingdi will register an increase in their population. This long-range phenomenon should be explained as a cascading effect, in which people from Cox's Bazar will begin moving to other south-eastern districts, prompting migration of local residents toward neighboring districts and ultimately gravitate toward Dhaka, Rajbari, and Narsingdi.

SLR Will Trigger Migration Patterns Involving the Whole Country
In Figure 5, we compare the results of the model by Davis et al. (2018), with our predictions, starting from the same initial populations in 2050 and accounting for the same values for the affected populations in the inundated districts (a total of 816,000 people over 21 districts). The two models anticipate a comparable net change in the population of all the districts (992,000 vs. 1,354,000) and are in agreement in the identification of the districts which will see their population reduce the most (Narayanganj, Shariatpur, and Munshiganj).
However, the predictions of the destination of the migrations are radically different. The model by Davis et al. (2018), predicts that SLR will cause a massive flux toward Dhaka, which will increase its population of more than 207,000 people, and northern and south-eastern districts will be unaffected by the internal migration. In contrast, our mathematical model predicts that migration will involve the whole country, including the northern districts, and that no district will register an increase in population above 67,000 people. Our model also suggests that Dhaka will see its population reduce by about 34,000 people. Differences between the two model predictions are further clarified in Table 1  With respect to the model by Davis et al. (2018), the proposed approach presents the following differences. First, the model by Davis et al. (2018) assumes that environmental change will cause the migration of every individual in an affected region, without considering unwillingness and inability to move by many people. In our model, we take as input the fraction of people that could migrate in an affected region and predict the ultimate migration. Second, the model by Davis et al. (2018) distributes migrants to neighboring regions, without accounting for either return migration or cascading effects, which are key factors in driving our model to its new equilibrium. Third, the model by Davis et al. (2018) is not formulated in terms of a continuity equation in which migration fluxes are combined to calculate the variations in the population throughout the network. As a result, it is tenable that the population distribution associated with the approach by Davis et al. (2018) does not constitute an equilibrium, from which there will be additional influx or outflux of people due to the immediate or long-term consequences of environmental migration.
The latter difference between our work and the approach by Davis et al. (2018) could be further developed by examining the transient dynamics of our model. In Figure 6, we show the population in the four districts of Dhaka, Narayanganj, Shariatpur, and Munshiganj as a funcion of time; these districts are selected based on Table 1, which defines them as the most critical districts for our model and the one by Davis et al. (2018). In all of these simulations, we use as initial conditions the predicted population at the equilibrium in the absence of inundations and, at the first time-step of the simulation, we modify the vector α to capture the simultaneous SLR of the 21 inundated districts. Results presented in Figure 6 confirm the intuition that the district of Dhaka might undergo an initial increase in its population, which would then reduce as time progresses. The same response is not visible for the other three districts, which are directly impacted by inundations and steadily register a decrease in their population.
These predictions depend on the modality that we used to force the system to a new equilibrium, whereby changing the spatio-temporal patterning of SLR-driven inundations would lead to different transient dynamics. To explore this aspect, we examined two alternative inundation patterns. Within the first pattern, one district was inundated at each time-step, taking 21 time-steps to inundate the whole set of districts (a total of 100 simulations was conducted to account for the randomness of the selection). Within the second pattern, the 21 districts were gradually inundated over 21 time-steps, by linearly increasing the inundation from 0 to the desired final value corresponding to the original analysis. In general, both of these two modifications delay and mitigate the migration process with respect to the original analysis. For example, in the district of Dhaka, we observe that a gradual inundation would delay the peak in the population of 20 timesteps and reduce its magnitude by 32.14%. Randomizing the inundated districts has an equivalent effect, although the extent to which the peak is delayed and reduced depends on the specific temporal ordering of the inundated districts. A slightly different scenario is recorded for the three other districts of Narayanganj, Shariatpur, and Munshiganj, where any of the two modifications is responsible for a diffusion of the migration over a wider time-window.

Summary and Concluding Remarks
Borrowing the words in Jasso (2003), "Every year and at every age humans move." They move to pursue a better quality of life, as the Romans put it, "ubi bene, ibi Patria" ("where one is well-off, there is one's country"), sometimes just to survive. Environmental change from droughts, desertification, floods, earthquakes, and wildfire is being increasingly recognized as a key driver for mass migration (Kaczan & Orgill-Meyer, 2019). SLR, in particular, is expected to affect the lives of more than 600 million people living close to the coast, who could see their livelihood at risk and decide to migrate (Hauer et al., 2019 Notes. Comparison between predictions by Davis et al. (2018) and our proposed modeling framework; the rank is based on the absolute value of the variation. An asterisk identifies a non-inundated district.

Table 1 Net Variations in the Population Due to Sea Level Rise in 2050 in Selected Districts
Here, we make a further step toward a universal model of human migration under environmental change. Our migration model extends the promising approach of Davis et al. (2018) along several directions to anticipate salient phenomena that have been extensively documented in the literature: unwillingness or inability to migrate, cascading effects on migration patterns, and return migration (Hauer et al., 2019). These extensions do not come at the expense of increased complexity, whereby the model contains a single additional parameter that is used to capture the baseline tendency of people to migrate, in the absence of environmental change. This is the only parameter that should be calibrated, before embarking on the predictions of migration patterns.
In principle, the modeling framework is applicable to a wide class of environmental migration problems. Similar to Davis et al. (2018), we demonstrate the framework on SLR in Bangladesh, using easily available geographic and demographic data to draw predictions on internal migrations. Our findings point at a more complex interplay between geography and demographics in shaping migratory fluxes due to environmental change than Davis et al. (2018). Highly populated districts should not always experience increases in population in response to environmental change occurring in lesser populated regions, but may also see their population reduce. Similarly, lesser populated regions could register further reduction or increase in their population, depending on where the location of environmental change. We propose that the predictions of our model should be preferred to those by Davis et al. (2018) in the study of SLR, whose impact will gradually affect the Bangladesh population, rather than drive immediately the system out of equilibrium.
The present model is not free of limitations, which should be explored in further studies.
1. During each time-step, we always assume that the hypotheses of the radiation model by Simini et al. (2012) hold true. In particular, we hypothesize that only economic factors drive the selection of the migration destination, thereby neglecting social ties which could play a critical role in shaping migration patterns as evidenced by Kartiki (2011), for rural Bangladesh. This limitation might be addressed by informing the construction of the weighted network of the radiation model based on social ties within communities. This modeling route could also assist in a more detailed quantification of the opportunities that are offered by large cities, against rural areas. Whether these improvements will result in substantially different predictions is presently unclear, whereby preliminary analysis of the heterogeneity in the baseline migration propensity seems to point at robustness of our model predictions against these unmodeled effects 2. In its present form, the model does not distinguish between recent and past migrants in each district, which may be crucial in determining the real opportunities available to migrants in each district. This DE LELLIS ET AL.
10.1029/2020EF001931 12 of 14 limitation might be addressed by extending the model beyond a first-order difference equation to capture sedimentation in the migration process. In this vein, one should explore the possibility of including memory in the migration rates 3. Introducing a memory in the migration rates may also help resolve the ambiguity underlying the notion of return migration. In its present form, the model accounts for return migration but does not offer a means to tease out their net effect. Specifically, from the flux between two districts, it is impossible to identify individuals who migrate for the first time versus individuals who are returning to their district of origin 4. The model does not distinguish between capability and vulnerability as two different causes for the decision to migrate (Kaczan & Orgill-Meyer, 2019). The complex interplay between these drivers might challenge the premises of the radiation model, which does not distinguish between individuals in their decision process. This limitation might be resolved through a further degree of refinement that partitions the affected population in classes 5. The model does not account for within-region migration, which could play an important role in migration patterns. For example, in Bangladesh, it has been proposed that a relevant fraction of displaced people could attempt short-range migration (Fussell et al., 2014). This limitation might be addressed by teasing out the portion of individuals who decide not to migrate from those who choose to migrate within the same district through an additional parameter to be identified from experiments 6. The model does not allow for a detailed treatment of external migration, in the form of crossing the borders of the country through ground transportation or airplane. Preliminary analyses on external migration to a virtual destination that is equally accessible from any district point at a secondary effect of the resilience of the country to localized shocks. Future research should explore whether this claim generalizes to more accurate representations of external migration Whether the simplicity of the modeling framework should be regarded as a limitation or a strength depends on the envisioned application of a mathematical model of environmental migration. The proposed approach requires the calibration of a single parameter to capture baseline migrations, which avoids overfitting the data, but comes at the cost of a limited temporal and spatial resolution of baseline migration throughout a country. Similarly, the prediction of environmental migration patterns is based on the number of individuals in any region who could become mobile due to environmental change. In principle, this number could account for economic, political, social, and demographic factors that could identify the most impacted individuals (Black, Adger, et al., 2011). However, working with the Bangladesh case study, we chose to estimate this number based on the extent to which each district will be inundated. This simplistic choice allows for a first assessment of the effect of environmental change, upon which one can determine the effect of specific, intervening factors.

Data Availability Statement
Data are available from the supplementary data published by Davis et al. (2018).