FeliX 2.0: An integrated model of climate, economy, environment, and society interactions

The Full of Economic-Environment Linkages and Integration dX/dt (FeliX) model is a System Dynamics-based Integrated Assessment Model (IAM), explicitly incorporating human behaviors and their dynamic interactions among global systems. This paper presents FeliX 2.0, describing the detailed framework and key interactions among nine integrated modules. FeliX 2.0 refined its original version in population dynamics, food and land use systems, and socioeconomic settings for poverty analysis. Robust calibration is applied to key variables against their historical data since 1950. Future projections of multiple variables up to 2100 demonstrate coherences between FeliX 2.0 and the IAMs used in IPCC assessments. Both outputs (the robust calibration results and future projections) underscore the efficacy of FeliX 2.0 in capturing complex interdependencies within global systems. FeliX 2.0 stands as an informative tool and offers insights into interactions within the human-Earth system and the analysis of complex economic-environmental-social challenges in short-and long-term future.


Introduction
Integrated Assessment Models (IAMs) are computational tools designed to simulate linkages in economic-environmental-social systems (Nordhaus, 1993;Schneider, 1997;Schwanitz, 2013;Van Beek et al., 2020).They aim to capture and represent in a simplified way to encompass a multitude of interconnected variables from diverse disciplines, such as greenhouse gas emissions, energy systems, land use, and economic growth.The greatest strength of IAMs lies in their ability to facilitate scenario planning, enabling decision-makers to explore various policy alternatives and understand the trade-offs involved in each option (IPCC, 2022a).By far, IAMs have been widely used for modeling and simulating global challenges such as climate change (IPCC, 2022a;2022b;Lontzek et al., 2015), sustainable development (Moallemi et al., 2022a;Pedercini et al., 2020), and energy transitions (Pye et al., 2015;Van Vuuren et al., 2016).However, most of existing IAMs tend to overlook feedback perspectives and nonlinear interactions among different real-world systems, potentially introducing inaccuracies into their predictions (Beck and Krueger, 2016;Schneider, 1997).For example, most IAMs don't endogenize feedback between economy and society or simplify human behavior as stylized demand assumptions driven by externally imposed macro-demographic and socioeconomic trends and/or models (Calvin and Bond-Lamberty, 2018;Van Vuuren et al., 2012).Uncertainties in model inputs and assumptions further complicate their outputs, and the extensive data requirements can be a hurdle in their applications.
A particular way to address IAMs' limitations in recognizing nonlinear interdependencies and feedback within real-world systems is applying System Dynamics (SD) modelling.SD is a well-established methodology with notable applications in understanding and managing systems that change over time (Meadows and Wright, 2008;Sterman, 2000).Widely recognized for its effectiveness, SD has been used in modeling feedback interactions, delayed responses, and non-linear behaviors.Its applications are not only in decision sciences but also within the broader context of the climate and sustainability (Moallemi et al., 2021).Illustrative examples of SD models include the World3 (Meadows et al., 1972), the Integrated Sustainable Development Goals (Millennium Institute, 2021), and the Earth4All (Dixson-Declève et al., 2022).These models simulate how various policies affect human wellbeing, societies, and ecosystems in the short and long term.Two unique features have been recognized, especially for climate-and sustainability-oriented SD models.One is the endogenous modeling of interactions between social, economic, and environmental systems in one integrated model.This feature is in line with the main purpose of IAMs, which seeks to recognize the interlinkages between economic-environmental-social systems.The other feature is the development of representations where relationships between observed outcomes and modeling assumptions are transparent and relatively easy to understand.Aggregate and descriptive SD models can provide valuable complements to insights derived from other Earth system models and detailed IAMs (e.g., partial or general equilibrium) that focus more on a detailed perspective on biophysical and socioeconomic systems but may lack the explicit modelling of feedback between systems.
The Full of Economic-Environment Linkages and Integration dX/dt (FeliX) model is a SD-based IAM that simulates complex interactions among global systems, including population, education, economy, energy, water, land, food, carbon cycle, climate, and biodiversity.The FeliX model is one of the very few models that explicitly incorporate human behavior in the human-Earth systems, considering comprehensive dynamic interactions of socio-economic and environmental sectors.It addresses a key limitation of conventional IAMs (i.e., the limited feedback and nonlinear interactions among systems) and covers the breadth of social, economic, and environmental aspects, and their interactions, in one integrated model.The FeliX model is a continuous model, operating on an annual time scale, and is designed to project global-scale future socioeconomic development and environmental conditions over the long term up to 2100.All the modules in the FeliX model have been calibrated to the historical data.This ensures the reliability of FeliX model in reproducing historical trends and producing plausible projections for the future.
This paper aims to present a systemic and comprehensive framework of the FeliX 2.0 model (shortened to FeliX 2.0 from here on) and describe key interactions among system elements included in FeliX 2.0.The remainder of this paper is organized as follows: Section 2 provides a brief overview of the previous research that used FeliX.Section 3 describes the model by elaborating on variables, functions, and interactions among different modules of FeliX 2.0.Section 3 presents the results about fitness of FeliX 2.0 outputs against historic data, and a comparative analysis of the model's future projections against those of other existing IAMs.Limitations, outlooks, and conclusions of FeliX 2.0 will be summarized in Section 4.

Overview of the previous research with the FeliX model
The initial version of FeliX model was presented by Obersteiner et al. (2012) and Rydzak et al. (2013), providing a conceptual description of the dynamic interactions and feedback loops of the model, along with mathematical formulations for all model variables.Based on the FeliX model, researchers have conducted scenario analyses on different global socio-economic and environmental problems.Thereafter, two main research directions, climate change mitigation and sustainable development, have been explored using FeliX.
Climate change mitigation.This direction aims to narrate emission pathways based on different mitigation measures.Walsh et al. (2015) analyzed emission pathways if microalgae was used as a feedstock for livestock and found that up to 2 billion hectares of land currently used for pasture and feed crops could be freed.In a subsequent study, Walsh et al. (2017) projected energy and land-use emissions mitigation pathways through 2100 and found that anthropogenic emissions need to peak before 2025 to maintain realistic pathways to meeting the COP21 emissions and warming targets.Thereafter, the model was extended to FeliX 2.0, focusing on population dynamics, education levels, food and land use systems including explicit fertilizer use and behavior-driven diet shifts, and socioeconomic settings for global poverty elimination.
Sustainable development.This direction explores how FeliX 2.0 can be employed in assessing the progresses of Sustainable Development Goals (SDGs).Eker et al. (2019) extended the FeliX model by developing a behavioral diet shift module in addition to the food, land use, and agricultural fertilizer consumption modules, and identified the main drivers of global diet change and its implications for the food system.Eker et al. (2019) showed that the social norm effect (for instance, the extent of vegetarianism in the population that accelerates a further switch to a vegetarian diet) and self-efficacy are the main drivers of widespread dietary changes.Moallemi et al. (2022b) implemented the global change scenarios in FeliX 2.0 to explore the impacts of model uncertainty and its structural complexity on the projection of sustainable development under these scenarios.Furthermore, Moallemi et al. (2022a) explored the drivers of progress of global SDGs and found that early planning for systems change into more sustainable pathways was important for accelerating the progress toward increasingly ambitious sustainable targets.Liu et al. (2023) further extended the FeliX 2.0 model with a global poverty module, to assess the effectiveness of different socioeconomic and environmental policies in eliminating poverty and quantified their impacts on the Earth system by 2050.More recently, an indicator for measuring human wellbeing, the global years of good life (YoGL) indicator, was estimated using the FeliX 2.0 model and found to increase by on average about twice the rate by 2100 (Kuhn et al., 2023).The model has been used in further applications to explore sustainable development pathways (Ruan et al., 2023), poverty alleviation strategies (Liu et al., 2023b), and food system transformation (Yang et al., 2023).

Structure and module descriptions of FeliX 2.0
The structure of FeliX 2.0 represents relations and interactions between nine real-world global systems: population (including education), economy, energy, water, land, food (including diet change), carbon cycle, climate, and biodiversity.The model outcomes are determined by many interacting feedback loops within and between these systems as shown in Fig. 1 and Table 1.
FeliX 2.0 consists of nice modules, that is, population and education, economy, energy, water, land use and fertilizer use, diet change, carbon cycle, climate, and biodiversity, to stimulate the ten real-world global systems.Key interactions among different modules are listed in Table 1.

Population and education
In contrast to many IAMs that take population projections as an exogenous input based on different Shared Socioeconomic Pathway (SSP) narratives (Kc et al., 2018;Riahi et al., 2017), the population module is endogenous and conceptualized based on dynamic mechanisms of population development and population ageing.Population development is governed by the population growth mechanism and the balancing effect of economic growth and education.The growth mechanism captures how a growing population leads to a higher economic output and how the economic growth increases the life expectancy at birth, reduces mortality, and hence further increases the population.The balancing effect mainly describes effects of increasing economic outputs and educational attainment decreasing fertility rates.The interaction of the feedback mechanism allows for simulation of population ageing.

Ageing structure
Population growth is based on an ageing chain and computes the male and female population size of 5-year age groups between 0 and 100+ years old (Fig. 2).The ageing chain represents the transition of newborns through the age cohorts as they age.It means that each age cohort (except the "0-4" cohort) has one inflow (maturation of the previous cohort) and two outflows (maturation to the next cohort and Matur gender,age is formulated as the division of population by the interval duration (Intvl Dur, 5 years).It assumes a homogenous distribution of population within the age group.

Birth rate and fertility
The birth rate (Birth gender , Eq. ( 2)) is the main factor affecting population dynamics, alongside the reproductive female population represented by gender and age-cohort segmentation.The birth rate is driven by education, and wealth represented by gross world product (GWP) per capita (GWP per Cap).Centre, 2020).This formula of logistic functions is extensively used in FeliX 2.0 to estimate impacts of a certain factor on interested variables.In consideration of the space limit in the main text, we refrain from providing a detailed elaboration of each logistic function from here on.Details about the logistic functions are available in the supplementary model documentation (Ye and Eker, 2024).
Tot Fert represents the number of births per woman at reproductive Finally, LE is converted to life expectancy for each gender and age group (LE gender,age ) with constant coefficients, which are estimated as the mean of historical values between 1990 and 2020.Associated data is obtained from the Global Burden of Disease dataset (Lopez and Murray, 1998).

Educational attainment and mean years of schooling
3.1.4.1.Educational attainment.Population by age cohort serves as input into the education module to compute the population of primary, secondary, and tertiary education graduates according to the enrollment rates and graduation rates (Fig. 1).Population with each educational attainment level, similar to the population chain, accounts for the ageing of people who graduate from each level and for transitions between the education levels.In detail, primary, secondary, and tertiary education graduates accumulate by each gender and 5-year age group corresponding to the respective education level.The assumptions regarding the age groups corresponding to education levels are as follows: children aged 5-9 enroll in primary education, with an average duration of 6 years.Subsequently, primary graduates, aged 10-14 or 15-19, may proceed to secondary education for another 6 years.Tertiary education enrollment is limited to the "15-19", "20-24", and "25-29" age groups, with an average duration of 5 years.Net change rate of primary education graduates (dPEG/ dt) for each gender-age group is formulated as: where Grad Prim gender,age is graduation rate from primary education, which is equal to enrollment rate to primary education (Enrl Prim gender ) yet after a delay of average duration of primary education.Enrl Prim gender is calculated by multiplying the population at the age of 5-9 with the primary education enrollment fraction for each gender (PEF gender ).PEF gender is defined as a reference value (Ref PEF gender , as the maximum possible enrollment fraction) multiplied by the impact of wealth (Imp GWP on Prim gender ).Imp GWP on Prim gender is formulated as a logistic function and calibrated based on GWP per Cap.
Matur Prim gender,age and Death Prim gender,age in Eq. ( 7) are maturation and death rates of primary education graduates, respectively, which is proportionally to the ratio between PEG gender,age and Pop gender,age .
The formulation of secondary and tertiary education graduates follows Eqs. ( 7) and ( 8), with differences in the age groups explained in the list of assumptions above.Enrollment to secondary (tertiary) education is assumed to be fraction of the primary (secondary) education graduates, implying that the previous education level is a prerequisite for enrollment.The effect of GWP on enrollment for each education level is calibrated to historical data from WCDE for the period 1950-2020 (Wittgenstein Centre, 2020).

Mean years of schooling (MYS).
Mean years of schooling serves as a crucial metric in delineating the effect of education on fertility rates and life expectancy.Mean years of schooling is formulated as the weighted average of the duration of each education level.The weights are determined by the total number of graduates for the respective education level relative to the population aged 15 and above.

Labor force
Labor force is distinguished as skilled (LF skilled ) and unskilled labor force (LF unskilled ).The skilled labor force is the sum of total population aged 15-64 with tertiary education and half of the population aged 15-64 with secondary education.The unskilled labor force is determined by the remaining population aged 15-64.

Economy
The economy module is based on the Cobb-Douglas production function, where Gross World Product (GWP), that is, global total GDP, is computed from labor input, total capital input from the energy and non-energy sectors, and total factor productivity from energy and nonenergy technologies.The Cobb-Douglas function is also incorporated the impacts of changes in ecosystems and climate change on the economic outputs.The assessment of global poverty rate is also included in the economy module.

Gross world product (GWP)
GWP is determined by the total reference economic output (REO).The total REO is the sum of the REO generated by the labor force by skill (REO skill ).Skilled and unskilled labor force are distinguished.REO skill are computed based on a Cobb-Douglas production function, depending on the technology (φ Tech skill ) and capital (φ K skill ) allocated to the labor force, and the size of the labor force.
where CEO skill is the capital elasticity output for skilled and unskilled labor force.REO Init

Climate change impacts on GWP
Impacts of climate change on economy (Imp CC on GWP) assumes that the temperature change has a direct impact on the economic output, following the commonly used 'damage function' forms presented by Nordhaus (2017), Dietz and Stern (2015), Burke et al. (2015), and a custom function defined in a flexible logistic form.The baseline damage assumption follows the estimates of Burke et al. (2015) which reach almost 75% of GDP loss at 5 • C of warming for the case where long-run (5 year) effects of temperature change are considered and countries give differentiated adaptation response.Details about the damage functions   (DFs) are available in (Eker et al., 2023) and Ye and Eker (2024).

Poverty
The global poverty rate (PR) is defined as the proportion of the population aged 15+ living below the international extreme poverty line (PL).The global PR is calculated as the sum of the poverty rates across the different population groups aged over 15 (PR gender,age ) weighted by their corresponding population shares.In addition, PL is defined as minimum consumption or income of a person per day and set as $2.15 per capita per day in 2017 PPP by The World Bank (2022  (11) where Real Income Param gender,age is determined by model calibration according to the data from The World Bank (2020).

Energy
The energy module includes total energy demand and energy supply from different fossil fuels (oil, gas, coal) and renewable energy (wind, solar, biomass).

Energy demand
Global energy demand is determined by the per-capita energy demand and the global population.The energy demand per capita (En Dem per Cap) is determined by the wealth of the population.

En Dem per Cap(t) = Imp GWP on En Dem(t) × En Dem per Cap Max
(12) Imp GWP on En Dem represents the impact of GWP on En Dem per Cap.It is modeled as a non-linear relationship-initially increasing at a low rate, then at an increased rate and eventually flattens out over the long run.Furthermore, Imp GWP on En Dem is calibrated according to the historical energy demand pattern obtained from Energy Institute (2023).En Dem per Cap Max is the maximal reference of energy demand per capita.

Energy production
Energy production comes from three fossil fuels (coal, oil, and gas) and three renewable (solar, wind, and biomass) sources, to meet the total energy demand.The production of any type of energy is determined by its identified resource and production capacity, and its market share.
3.3.2.1.Energy production from fossil fuels.The production of fossil fuels is modeled as a supply chain structure, adapted from Sterman and Richardson (1985) and Davidsen et al. (1990).This section uses oil as an example to explain main components of the energy production structure for fossil fuels.In general, oil is produced from identified oil resources, which grows based on oil exploration through existing technologies from undiscovered oil resources.
Oil exploration (Expl oil ) is counted as the lower value of a potential exploration rate of oil (Expl Pot oil ) and a desired exploration rate of oil (Expl Des oil ).Expl Pot oil represents the exploration possible due to effective investments in oil exploration (Eff Expl Inv oil ) and productivity of available technologies (Prodvty Expl Inv oil ) which encapsulates the resource depletion effect.
Expl Des oil represents the amount of oil resources that would ensure continued production at the current rate.Unit cost of oil exploration is then estimated depending on the remaining undiscovered oil resources and advances in exploration technologies.The cost of oil exploration in turn determines the desired investment in oil exploration.
Similarly, the oil production is counted as the lower value of either total oil demand or a potential oil production rate.The total oil demand is calculated as the market share of oil in total energy demand.The potential oil production rate equals to the production due to effective investments in oil production and available resources.Similar to the costs of exploration, the unit cost of oil production is then dependent on the productivity of investment in oil production which encapsulates the resource availability effect.It is important to note that production costs also include the policy-set carbon price.The cost of both oil exploration and production determines the oil price, which affects the market share of oil demand in total energy demand (see section 2.3.3Market shares of energy sources in total energy demand).

Energy production from renewable sources.
There are several main factors that affect the production of renewable energy from wind, solar, and biomass sources.They are resource characteristics (e.g., average sun radiation for solar energy), efficiency and aging of technology, land availability, installation cost, investment in installation and efficiency, and the impact of technological learning on unit costs of installation.This section uses solar energy as an example to explain important dynamics in the structure underlying renewable energy production.
Firstly, total solar demand equals to the market share of solar energy multiplied by the total energy demand.Secondly, the solar production capacity is calculated to meet total solar demand but constrained by the potential solar energy production.The potential solar energy production depends on the available solar resource and technological constraints of average intensity of the resource, annual availability of the resource considering weather impacts, changing efficiency of solar radiation conversion at the current state of technical developments, and the installed capacity to convert solar radiation into energy.The installed solar capacity changes according to the new installation of solar capacity, which is influenced by the productivity of solar investments and the adjustment of existing solar infrastructure.Lastly, the unit cost of solar energy is estimated as a sum of the unit costs of installation and production.

Market shares of energy sources in total energy demand
The total energy demand is allocated to different energy sources based on their respective market shares (MS).MS is calculated based on indicative market shares (MS Ind) of energy sources which are defined to represent accumulated market shares of different energy sources over time.The price of each energy source is averaged over a period, and all the prices result in the average energy price.The price competitive factor is calculated for all energy sources as a ratio of the averaged price of each energy source and the average energy price.The price competitive factor is further adjusted by the elasticity of that energy source price to demand.

Water
The water module mainly quantifies global average water scarcity, that is, the balance between water availability and use.Water availability is a function of available water resources, water recovery of used water in economic sectors, and a drought rate for the impact of climate change.Total water use considers water withdrawal by three economic sectors (agriculture, industry, and domestic) to meet their water demands.

Water supply and water availability
Available water resources consist of water supply rate of fresh and non-conventional water, water recovery rate of used water resources by the three economic sectors, and drought out rate representing extreme drought events.Water supply is estimated by multiplying desired water supply rate and water supply fulfillment rate (describing water supply fulfillment as a relation of water supply and demand).Desired water supply rate ( Wat Sup Des) is to cover water consumption (Wat Cons, i. e., water use decreased by recovery of used water resources of three economic sectors) and available water resources adjustment (Wat Avail Adj, taking into account dynamics of total water demand and water safety stock coverage).

Wat Sup Des(t) = MAX(0, Wat Cons(t)) + Wat Avail Adj(t) (15)
With growing water demand, the supply fulfillment might be impaired which relates to infrastructure design and its operating limits.This relationship is modeled as a logistic function dependent on the amount of water that can be reliably provided on annual basis for agricultural, industrial, and domestic sectors due to resources and infrastructure availability.Associated data are obtained from 2030 Water Resources Group (2009).

Water demand, water withdrawal, and water use
Water is demanded for economic production by agriculture, industrial, and domestic sectors.Agricultural water demand (Wat Dem Agri) is dependent on two ways of land watering-irrigation and rainfed.
where Wat Int Irr and Wat Int Rfed represent average agricultural water demand per square (i.e., water intensity, Wat Int) agricultural land areas of irrigation (Area Irr) or rainfed (Area Rfed), respectively.Area Irr is calculated by multiplying total agricultural land areas with the percentage of irrigated land areas, which is estimated by considering the impact of wealth (represented by GWP per Cap) and their upper and lower limits.This approach is also applied to estimate Wat Int Irr, Wat Int Rfed, and total water demand for industrial sector based on their own upper and lower limits, respectively.Domestic water demand is quantified based on total population and average domestic water demand per capita.Domestic water demand per capita is estimated by considering the impact of wealth (represented by GWP per Cap) and their upper and lower limits.
Agricultural, industrial, and domestic water demand drives water resources withdrawal from available water resources.Similar to water supply, agricultural, industrial, and domestic water withdrawal rates slow down when their demands approaching max water withdrawal rate.This process is represented by water withdrawal fulfillment rate (Fulf Rat Withd).Particularly for agricultural water withdrawal (Wat Withd Agri), impacts from extreme drought (Drght Rat) is also considered to simulate the decrease in Wat Withd Agri and available water resources due to drought.
Agriculture, industrial, and domestic water withdrawal rates accumulate as total used water resources.As mentioned in the section Water supply and water availability, a fraction of total used water resources can be recovered and supplied as available water resources.

Land use and fertilizer use
The global land use and land use change dynamics are represented based on four main categories of land use as defined by FAO (FAOSTAT, 2020a): agricultural, forest, urban/industrial, and the other land that does not fall into any of the first three categories.Land use change refers to the conversion of land among these four categories.Considering the historical and expected dynamics, a bi-directional conversion between the agricultural land and forest land and between the agricultural land and other land is assumed.The rest of land use changes are considered one directional conversion, that is, from agricultural land to urban and industrial land, from forest land to urban and industrial land, and from other land to forest land (Fig. 4).

Land conversion
The main underlying driver of land use change is the food system, in addition to bioenergy and forest management practices (Fig. 4).Each land category is modeled as a stock.Apart from urban and industrial land, which is set to increase only, all other types of land can gain more space or lose it due to land conversion.For instance, abandoning any human activities in a part of agriculture land and allowing certain duration of time the agriculture land will be covered by grass, later by shrubs, yet later by trees and eventually it will be classified as forest.The easier it is for the natural process to transform the area the transformation time is set to a lower value.In case of forest-agriculture land conversion, agriculture-other land conversion, and other-forest land conversion, there are additional forces apart from natural conversion processes that drive and increase the rate of expansion or shrinking of the forest, other and agriculture land areas.A constant ratio is set for other land and forest land areas as protected land, respectively.In case of forest land, the protected area increases successively.The protected areas are excluded from any transformation processes.The agricultural land requirements and its estimates will be especially elaborated in the next section due to the key link with diet change module (Fig. 1).

Agricultural land requirements
Agricultural land requirement increases due to the growing population and income levels which lead to more food demand (Fig. 5).Agricultural land is further divided into three sub-categories as arable land, permanent cropland, and permanent meadows and pastures.Land use shifts due to food demand are considered based on plant-and animal-based food that have their own land allocation and yields.Eight categories of food products are specified and sorted into plant-based and animal-based food.Plant-based food include pulses, grains, vegetables and fruits, and other crops such as sugar and oil crops.Animal-based food are pasture-based meat, crop-based meat, dairy, and eggs.This categorization of especially animal products were motivated by the land use shares: The pasture-based meat category includes the red meat sources such as beef, sheep and goat, since 96% of the global average land use per unit of beef production is attributed to pasture land, that is, 12 of 12.5 ha per million kcal (Ranganathan et al., 2016).As for poultry and pork, the average land use footprint is less than 2 ha, while a large portion of this is on cropland (Ranganathan et al., 2016), since grains provide the 71% of the total feed demand (FAOSTAT, 2020a).
For each plant-based food, the arable land needed (ALN) is formulated according to the desired production rate (Prod Des) and the expected crop yield (Yield Exp), which is a 1-year time-averaged value of the crop yield over time.

ALN food (t) =
Prod Des food (t) Yield Exp food (t) , where food ∈ [pulses, grains, vegetables and fruits, other crops] Prod Des food of each plant crop is the sum of the demand for that crop to be used as food, livestock feed, and in other sectors.The estimate procedures of crops to be used as food and livestock feed are described in section 2.6 Diet Change.Crop to be used in other sectors is quantified as a constant ratio of total plant-based food demand.The actual production of each plant crop is the lower value of Prod Des food and the production volume calculated as crop yield multiplied by the area harvested.Crop yield is formulated with respect to a reference yield values in 2016 with the dynamic multiplicative impact of fertilizer application, water withdrawal rate and technology change.The total agricultural land requirements are estimated by dividing ALN food by the annual share of arable land in the total agricultural land.The rest of agricultural land is proportionally distributed as permanent cropland and permanent meadows and pastures according to their annual share in agricultural land.
The discrepancy between agricultural land requirement and available agricultural land pushes agricultural production systems either towards more fertilizer use to increase the crop yields and reduce the land requirement, or to deforestation and other land requisition to expand the agricultural land.Eventually, food supply is dependent on the land available, hence harvested, for each food category, and the crop (or meat) yields that depend on fertilizer use, water availability and other agricultural management practices.

Fertilizer use
Nitrogen and phosphorus fertilizer use in agriculture, from commercial sources or livestock manure, is explicitly modeled.Both nitrogen and phosphorus use, as well as their impacts, are modeled in a similar framework that focuses on chemical inflows to and outflows from agricultural land.This section uses phosphorus (P) as an example to explain the framework.
Global P balance with major inflows and outflows (Fig. 6) is based on an input and output framework (MacDonald et al., 2011) and a former modeling study about the biochemical processes of phosphorus (Dumas et al., 2011).In general, the global P residual in soil grows due to the application of commercial P fertilizers and the application of manure as a fertilizer, and decreases due to the uptake by crops and the loss of P via erosion, leaching to freshwater systems and surface runoff.The residual P and flows are reported in terms of the volume of elemental P, except the commercial fertilizer production and application, which are  reported as the volume of P 2 O 5 content.The application of commercial P fertilizers is modeled as the minimum of P demand from agriculture and P supply for agriculture.P demand from agriculture is formulated based on a reference application rate, taking into account the effect of income levels and land scarcity.Both effects are formulated as logistic functions based on GWP and the ratio between global agricultural land demand and the available agricultural land, respectively.Supply for agriculture is a fraction of the total P via P rock extraction.This fraction is set as 0.9, since the agricultural use has been on average 90% of the total P 2 O 5 produced between 2002and 2015(FAO, 2018)).P uptake by crops are estimated as a weighted average of the P contents in individual food category, using the data from USDA Food Composition Databases (USDA, 2020).We assume a nonlinear relation between P uptake and crop yields, based on the agronomic nonlinear dynamics of crop growth (Yin and Struik, 2010) that result in sigmoid growth.

Diet change
Food demand is quantified based on the total caloric demand for eight main food categories.Total caloric demand is based on dietary choices of different population segments.Population segments of dietary choices are represented by the followers of the meat-based and vegetarian diets for each age cohort and gender (Eker et al., 2019).The shifts between these two dietary choices depend on income (represented by GWP per Cap) and social and behavioral factors such as climate and health risk perception, self-efficacy and social norms that underly pro-environmental behavior.
Followers of meat-based and vegetarian diets are assumed to eat a standard mix of eight food categories.Typical United States' diets (Pimentel and Pimentel, 2003), a reference world diet and the global average supply statistics, were used as benchmarks to understand how meat-based and vegetarian diets differ globally.Based on such proportional differences of meat-based and vegetarian diets, the reference meat-based and vegetarian diets were formed by decomposing the world's average diet according to the population fraction of the two groups.

Caloric demand
Total annual caloric demand for plant-based food (Cal Dem) is the sum of total plant-based caloric intake by the vegetarian population (Pop Veg) and meat-eating population (Pop Meat).
Cal Dem food (t) = Intk Veg food (t)×Pop Veg(t)+Intk Meat food (t)×Pop Meat(t), wherefood ∈ [pulses, grains, vegetablesand fruits, other crops] ( where Intk Veg plant and Intk Meat plant are the annual per-capita intake of plant-based calories in the vegetarian diet and meat-based diet, respectively.Both are assumed as a constant fraction of the annual total caloric intake per capita.Annual total caloric intake per capita is formulated as the multiplication of a reference intake per capita and the effect of income.The latter is a logistic function based on GWP per Cap.Total annual caloric demand for animal-based food is simplified as where Intk Meat meat is the annual per-capita intake of meat-based calories in the meat-based diet.

Food demand and food production
Food demand is calculated as Cal Dem food divided by the caloric value of associated food category.In addition, plant-based food categories are also be used as livestock feed.The amount of plant-based food used for feed depends on crop-based meat demand, feed share of four plant-based food, and unit-feed requirement for livestock.Food pro-duction is eventually steered by waste-adjusted food demand, which is calculated as food demand divided by a waste fraction of supply.The waste fraction is set as 0.30 for grains, 45% for pulses, vegetables and fruits the waste fraction is 45%, and 20% for other food categories.
The production of animal-based food products (pasture-based meat, crop-based meat, dairy, and eggs) follows a similar approach based on grassland availability, food demand and land yield.Dairy production is aligned with meat production, and the production of eggs is formulated with respect to crop-based meat production, which depends on the feed fraction of crop production.

Carbon cycle
CO 2 emissions are calculated based on representations of carbon emissions from the energy and land-use change sectors.These emissions accumulate in the atmosphere until they are absorbed into the biosphere, pedosphere or oceans based on C-ROADS (Sterman et al., 2012(Sterman et al., , 2013)).

Carbon emissions
Total carbon emission rates from the energy sector include carbon emission rates from fossil and renewable sources.Carbon emissions of a specified energy (C Emis En energy ) are based on the carbon intensity of energy production (C Int energy ), which is calibrated to historical emissions within the uncertainty ranges of the unit emissions of energy

Carbon stock
Stocks of carbon are considered in the atmosphere, biosphere, mixed ocean layer, and four deep ocean layers (Fig. 7).Carbon in the atmosphere (C Stk Atm) accumulates through total carbon emissions formulated above.Carbon in biosphere is captured in biomass (C Stk Biom) and soil (C Stk Soil).C Stk Biom includes carbon stock in leaves, branches, stems and roots, whereas C Stk Soil includes carbon stock in litter and humus.As the concentration of C Stk Atm rises, it forces increase of the uptake by ocean and biosphere.Carbon flux from the atmosphere to the biomass is modeled according to the formula in Wullschleger et al. (1995) and grows logarithmically as the concentration of C Stk Atm increases.The residence of C Stk Biom depends on average lifespan.The outflow of C Stk Biom is partitioned between carbon flux from the biomass to the atmosphere and to the humus according to humification fraction.The outflow from C Stk Soil is equal to its content divided by its average lifespan in the humus.The flux between C Stk Atm and carbon stock in the mixed ocean layer (C Stk Ocn) adjusts to an equilibrium that considers buffer factor, a measure of the resistance to atmospheric carbon dioxide being absorbed by the ocean surface layer.The buffer factor itself rises with the atmospheric concentration which decreases ocean absorption capacity.Deep ocean diffusion fluxes are modeled as a simple eddy-diffusion structure.

Climate
The climate module is based on the C-ROADS model (Sterman et al., 2012), which in turn refers to FREE model (Fiddaman, 2002) and DICE model (Nordhaus, 1992(Nordhaus, , 1994)).The Earth's radiation budget is constrained to the temperature change due to carbon dioxide (CO 2 ), methane (CH 4 ), nitrous oxide (N 2 O), halocarbons and other forcings (e. g., aerosols, O 3 , etc.) CO 2 emissions are endogenous variables as described in the section Carbon cycle.The rest of forcings are generated exogenously based on historical data and the future projections from Representative Concentration Pathways (RCPs) (Byers et al., 2022).The temperature change is governed by radiative forcings, feedback cooling due to outbound longwave radiation, and heat transfer from the atmosphere and upper ocean to deep ocean layers.

Radiative forcing
Total radiative forcing consists of CO 2 radiative forcing due to increasing concentration of carbon in atmosphere, and other forcings which includes variables representing forcings from CH 4 , N 2 O, halocarbons and other gases and aerosols.

Feedback cooling
Feedback cooling due to outbound longwave radiation governs feedback mechanism of the atmosphere and the upper ocean.The rate of cooling is determined by the climate sensitivity-a metric used to characterize the response of the global climate system to a given forcing.
It is broadly defined as the equilibrium global mean surface temperature change following a doubling of atmospheric CO 2 concentration.

Heat transfer
The heat transfer into deeper layers of the ocean is modeled as a function of the eddy diffusion, which controls the movement of carbon through the deep ocean.Four different layers of deep ocean are considered.

Biodiversity
The biodiversity module is a simple structure representing population carrying capacity.The population is represented by global and encompassing all biomes mean species abundance (MSA).MSA is increased by species regeneration rate and decreased by species extinction rate.In addition, MSA approaching species carrying capacity (Spec Capa) limits species regeneration and intensifies species extinction.This process is quantified as logistic functions and based on the ratio between MSA and Spec Capa, and the regeneration factor and the extinction factor, respectively.

Species carrying capacity
Spec Capa is calculated based on reference species carrying capacity (Spec Capa Max), representing maximum sustained population size, and influencing factors related to fertilizers consumption (Imp Fertz on Biodiv), biomass production for energy purposes (Imp Biom on Biodiv), climate damage (Imp CC on Biodiv) and land use change (Imp Land on Biodiv).Spec Capa(t) = Spec Capa Max × Imp Fertz on Biodiv(t) ×Imp Biom on Biodiv(t) × Imp CC on Biodiv(t) × Imp Land on Biodiv(t) (23) The four influencing factors are estimated by logistic functions.In detail, Imp Fertz on Biodiv is a logistic function of fertilizer consumption (including nitrogen, phosphate, and potash fertilizers).Imp Biom on Biodiv consists of impacts from forest agriculture biomass production, which are logistic functions of related land areas.Imp CC on Biodiv is adopted from climate impact on economy, as a logistic function of temperature change from preindustrial.Imp Land on Biodiv takes into consideration changes of agricultural land, forest land, and other land compared to their initial areas.

Calibration to the historic data
Simulation trajectories of key variables from different modules in FeliX 2.0 demonstrate a strong alignment with their respective historic data (Fig. 8).Fifteen key variables are selected from different modules based on the best availability of their historic data and their interdependencies with other variables in the whole FeliX 2.0.The model is calibrated for the period 1900-2022, mostly 1960-2022 depending on the data availability, and then projects to the future.The historic data for these variables are sourced from reputable repositories, including Wittgenstein Centre Human Capital Data Explorer (Wittgenstein Centre, 2020) for population, total fertility rate, education graduates; the World Bank Data (The World Bank, 2020) for GWP, GWP per capita, and global poverty rate; International Energy Agency (IEA, 2020) for energy demand; FAOSTAT (FAOSTAT, 2020b) for agricultural land, and total daily calorie supply per capita; IPCC (2014) for total C emission from energy sector; NASA GISS (2023) for temperature change from preindustrial; and European Environment Agency (EEA, 2020) for atmospheric concentrations of CO 2 .
Q. Ye et al.The R square (R 2 ) values between trajectories simulated by FeliX 2.0 and historic data of these variables are more than 0.93 except of agricultural land (R 2 = 0.81) and temperature changes from preindustrial time (R 2 = 0.84).For the agricultural land, the low R 2 value can be attributed to the assumption that land use change, i.e. deforestation, is captured at higher aggregation level between forest and agricultural land, and a constant fraction of the total agricultural land is allocated to cropland and grassland.As for temperature changes from preindustrial time, the historic data per se were modeled such as observed with ±10% range of uncertainty according to (EEA, 2020; IPCC 2021).The best fitness is observed for the total population, total secondary graduates, total tertiary graduates, GWP, GWP per capita, and atmosphere CO 2 concentration.Their R 2 values are all more than 0.99.The best fitness of these variables results from relatively more available data.As such, given the good reproduction of historical trends in FeliX 2.0, the model is able to capture the dynamic connections and interactions of the complex socio-economic and environmental system for simulating future development.
The comparison shows that FeliX 2.0's projections fall within the ranges of results from selected IAMs, except for the total radiative forcing and GWP from FeliX 2.0 without DFs.Future total radiative forcing in both FeliX 2.0 is lower than those in other models.This can be attributed to our choice to use RCP4.5 projections for the non-CO 2 GHGs in the baseline scenario of FeliX 2.0, in line with the current policies, while the SSP2 baseline corresponds approximately to RCP7.0.When considering DFs in economic outputs, future GWP would decrease by more than 50%.Whereas, if no DFs are considered in FeliX 2.0 either, future GWP would still be consistent with the SSP inputs to the other models (Dellink et al., 2017) which do not take climate impacts on economy into account.Total CO 2 emissions from the energy sector before 2080 in FeliX 2.0 are projected to be higher than those in other IAMs.The exclusion of carbon removal technologies such as carbon capture and storage under SSP2 is the main reason for the highest CO 2 emissions in FeliX 2.0.Coal production is relatively low while wind energy production is relatively high in FeliX 2.0 compared to the production in other selected models.It is because the market share of coal would shrink to 24% in 2100 while the share of coal in energy supply in GCAM would be 30-50%.In addition, the market share of wind would expand to 12% in 2100 in FeliX 2.0.FeliX 2.0 results in relatively higher market shares of renewable energy (especially wind), resonating with the finding in Jaxa-Rozen and Trutnevyte (2021) about underestimation of the fast development of renewables by most IAMs.Projection results of crop production and forest land in FeliX 2.0s are both in line with results from other IAMs.

Limitations and outlooks
As a model developed for the purpose of capturing complex biophysical and socioeconomic mechanisms at an aggregate level, FeliX 2.0 offers a comprehensive framework that addresses the limitations of conventional IAMs, notably the inattention to feedbacks and nonlinear interactions.The applications of FeliX 2.0 (see Section 2) in scenario analyses, emissions pathways, and sustainability assessments have proved its reliability as a useful tool for understanding global system dynamics and human well-being.However, FeliX 2.0 is also subject to several limitations.The first limitation arises from its global scale, preventing fine-scale analysis of socioeconomic-environmental dynamics.This hampers the effectiveness of FeliX 2.0 in devising targeted policies for localized challenges (e.g., water scarcity and deforestation) and hindering comprehensive solutions tailored to specific geographic contexts, given high heterogeneity across regions and nations.In addition, it also prevents FeliX 2.0 from adequately capturing the diverse contributions of specific regions, to specific global problems such as climate change.To address this limitation, it requires to add a regional dimension in associated modules in FeliX 2.0 such as population, economy, etc.In addition, each regionalized module should follow the global structure and be calibrated by historical data of each region, to make sure the consistency and robustness of generated results.
The second limitation is the relatively macro-modeling approach of different sectors in FeliX 2.0 that does not involve a high level of technoeconomic detail.This reduces FeliX 2.0's ability to adequately capture the diverse contributions of specific sectors to specific global problems.For example, the end use sectors leading to the primary energy consumption are not modeled endogenously in FeliX 2.0.The transportation sector, for instance, accounted for around 30% and 20%, respectively, of global total energy consumption and CO 2 emissions in 2020 (IEA, 2022).Hence, it is also critical for climate mitigation through transitioning to electric vehicles, investing in public transportation, implementing sustainable traffic planning, and decarbonization in shipping and aviation (Jaramillo et al., 2022).To address the second limitation, more key sectors and factors of human-earth system can be added into FeliX 2.0 while focusing on the feedbacks rather than technoeconomic detail, as the distinguishing feature of FeliX.Lastly, FeliX 2.0 is a publicly available model, yet it requires a licensed software (Vensim DSS) for full functionality.A fully open access version of FeliX 2.0 can enhance its accessibility, transparency, open review, and collaboration, that are critical for robust and inclusive research.

Conclusions
The development and utilization of FeliX 2.0 have enabled understanding of global socioeconomic-environmental dynamics and their multilateral interactions.The robust calibration results of key variables within FeliX 2.0 modules against historic data highlights the model's efficacy in capturing the complex interdependencies of global systems since the beginning of the 20th century as they have been observed.Moreover, future projections from FeliX 2.0 exhibit coherence across multiple variables, including emissions, energy sector developments, crop production and forest land, compared to other existing IAMs.These insights highlight the strengths of FeliX 2.0 while also signaling avenues for refinement to enhance its usefulness and policy relevance, offering a promising trajectory for robust long-term scenario simulations of global dynamics.The diverse applications of FeliX 2.0 in areas such as climate change, sustainable development, and global poverty elimination have provided evidence of its efficacy as a powerful tool for short-and longterm analyses of interested economic-environmental-social systems and existing global problems.However, FeliX 2.0 also have limitations that require future enhancements.Examples are its global scale and limited sector coverages.In summary, FeliX 2.0 stands as an informative tool for both SD and IAMs, and offers insights into the human-Earth system feedbacks.Overcoming its limitations will boost more potential of FeliX 2.0, enabling future research to design targeted policies and comprehensive solutions for a sustainable future.

Fig. 1 .
Fig. 1.Overview of FeliX 2.0.We include 11 components in the figure to better illustrate interactions within and among the modules.As such, education shows as a separate component from population, and fertilizer use shows separate from land use.Source: Moallemi et al. (2022).

Fig. 2 .
Fig. 2. Inflows and outflows to different population cohorts in the population module.Population cohorts also distinguish gender into male and female.According to the System Dynamics methodology, boxes represent stocks/accumulations of elements while the arrows represent in-flows and out-flows indicating the direction of movement of quantities between stocks. dPEG

Fig. 3 .
Fig. 3.The lognormal probability density function of income.The shaded area is equal to the poverty rate.
Q.Ye et al.

Fig. 7 .
Fig. 7. Carbon stocks and associated flux among the stocks.Carbon accumulation in deep ocean is taken into account at four distinct ocean layers, yet not visualized for simplicity.

Fig. 8 .
Fig. 8. Calibration results of key variables in FeliX 2.0.In the labels of ordinate, B stands for billions, Deg stands for degree, Dmnl stands for dimensionless, ha stands for hectare, Kcal stands for kilocalorie, Mtoe stands for million tonnes of oil equivalent, ppm stands for parts per million, T stands for trillions, and TonC stands for tonne of carbon.

Fig. 9 .
Fig. 9. Future projections of key variables under the shared socioeconomic pathway 2 (SSP2) in FeliX 2.0 and five existing IAMs.FeliX 2.0 Ref represents outputs from FeliX 2.0 without considering the climate impacts on economy, that is, the damage function (DF), while FeliX 2.0 with DF represents outputs from FeliX 2.0 with climate impacts.SSP2 projections by other models are obtained from the AR6 Scenario Database (Byers et al., 2022).

and maturation rate (Matur gender,age ) by gender and age.
Q.Ye et al.  mortality).Population of each gender and age interval (Pop gender,age , Eq. (1)) accumulates by three flows, determining the net rate of change (dPop/dt).These three flows are birth rate (Birth gender,age ), death rate (Death gender,age ),

which is the number of births per woman in a particular age group during a 5- year period. Fert age is formulated as a function of total fertility (Tot Fert).
where g gender denotes the birth gender fraction, i.e., gender fraction of registered newborns.Fert age refers to the age-specific fertility rate, ).
Q.Ye et al.

Table 1
Main interactions within and among modules in FeliX 2.0.The interactions are direct interactions, always from the row modules to the column modules.

Land use, fertilizer use and food production
Mortality rate and life expectancy at birth Death rate (Death gender,age ) refers to the number of people who pass away in each gender and age group per year.It is formulated as a fraction of the population of each group.Death gender,age (t) = Pop gender,age (t) × Mort gender,age (t) (5) where Mort gender,age is the mortality fraction by gender and age cohort.
Bressler et al. (2021)culated as a logistic function of global average life expectancy at birth (LE).LE is derived through a multiplicative function considering impacts of wealth, education, total food supply, and climate change on a reference value of LE-the observed life expectancy at birth in 2000.The impacts on life expectancy stemming from wealth, education, and total food supply are estimated by logistic functions based on GWP per Cap (The World Bank, 2020), mean years of schooling (Wittgenstein Centre, 2020), and total daily calorie supply per capita (FAO-STAT, 2020a), respectively.The last impact, climate change's effect on LE (Imp CC on LE), is estimated depending on future temperatures.Global temperature rise is the main driver of climate mortality-the death rate directly attributed to the impacts of climate change such as extreme weather events, rising temperatures, and related environmental changes (IPCC 2022b).Meanwhile, education is expected to mediate climate mortality.To quantify Imp CC on LE, climate change impacts on mortality (Imp CC on Mort) is estimated firstly, by taking into account both temperature and education and calibrated to the estimates ofBressler et al. (2021).Imp CC on LE is inversely proportional to Imp CC on Mort.

Table 1
(continued ) Tech Eng is endogenously determined by investments in the energy module while Tech Oth follows an exogenous trend and data.The labor force (LF gender,age,skill ) is the corresponding labor force by different skills multiplied by the labor force participation rates for the respective groups (by gender and age cohort), which is set to be 34-78% depending on genders and age cohorts following The World Bank (2020).
1900,skill , φ Tech skill , φ K skill , and CEO skill are exogenous parameters determined by model calibration based on historical data of GWP and GWP per Cap from The World Bank (2020).K and K Init are annual capital stock and initial capital stock in 1900, respectively.Technology-related factor productivity (Tech) is distinguished into energy technology (Tech Eng) and all other technology (Tech Oth).Particularly, Enrl Sec gender,age (t) − Matur Prim gender,age (t) − Death Prim gender,age (t); if age in {˝10 − 14˝, ˝15 − 19˝} Matur Prim gender,age− 1 (t) − Matur Prim gender,age (t) − Death Prim gender,age (t); if age > ˝20 − 24˝( Q.Ye et al. Income per Cap gender,age is related to the GWP per capita.Income per Cap gender,age (t) = GWP per Cap(t) × Real Income Param gender,age ). Income per capita by gender and age cohort (Income per Cap gender,age ) is used to calculate PR gender,age(Fig.3,Fosu,2010;Hughes, 2015).Income per Cap gender,age is assumed to follow a log-normal distribution characterized by the mean μ gender,age and standard deviation σ gender,age of the normal distribution function of ln ( Income per Cap gender,age ) , based on previous research (Fosu, 2010; Lakner et al., 2022; Liu et al., 2023a).PR gender,age (t) ( Income per Cap gender,age ≤ PL ) = ∅ ( ln(PL) − μ gender,age (t) σ gender,age (t) )(10)The value of a standard normal cumulative distribution function can be obtained by looking up the standard normal distribution table.
Eff MS Price energy (t) − MS Ind energy (t) Time Adjwhere MS Ref energy is the reference market share of a specific energy source.Time Adj is set as 10 years for all energy sources.Eff MS Price depends on the competitive price factors and prices of energy sources.
production(IPCC, 2014).C Emis En energy(t) = C Int energy × Prod energy (t)(21)where Prod represents the annual production of different types of energy.Carbon emissions of fossil fuels also include the effect of carbon capture and storage technology.Carbon emissions from land-use change (C Emis Land) includes emissions from agricultural land change (C Emis Agri), and forest land change due to deforestation and forest conversion to managed forests and plantations (C Emis Frst).Int Agri and C Int Frst represent carbon intensity of agricultural land use and forest land use, respectively.Land Chg Agri and Land Chg Frst represent ratio of agricultural and forest land area change compared to its initial area in year 1900.