The projections were created using the Proudman Oceanographic Laboratory Coastal Ocean Modelling System (POLCOMS; Holt and James, 2001) coupled to the European Regional Seas Ecosystem Model (ERSEM; Butenschön et al., 2016). Together, these simulate the movement of water, energy and dissolved and suspended material through the sea and the cycling of nutrients and carbon through the marine ecosystem.

POLCOMS is a three-dimensional model of physical processes, suitable for modelling both deep and shallow water and areas with steep bathymetry. It is a free surface baroclinic model, with varying water depth, and includes the effects of tides but not waves. Forty depth levels were used at each point, regardless of total water depth, distributed more closely in the upper parts of the water column than at depth (a modified sigma system). The minimum water depth was 10 m. When run at the 0.1° resolution used here, POLCOMS is able to explicitly model processes of riverine freshwater inputs, tidal mixing and coastal upwelling that are parameterized in CMIP5 global models (Holt et al., 2009; Drenkard et al., 2021). These processes can influence circulation patterns and nutrient distributions in near-shore waters, with potential effects on water temperature and productivity which could affect resources such as corals, aquaculture sites and production at higher trophic levels. Given the widely varying water depths and very steep shelf edges within the region ( Figure 1 ) it is important to use a model with adequate vertical resolution – large parts of the shelf area are less than 100 m deep.

ERSEM models the transfer of carbon, nitrogen, phosphorus and silicate through the lower trophic levels of the marine ecosystem (Butenschön et al., 2016). It is one of the more complex models of its type and is well suited to modelling coastal and shelf-sea environments. It has four phytoplankton functional types, three zooplankton types and bacteria, as well as six categories of dissolved and particulate organic matter. The carbonate system is included, enabling changes in pH to be modelled. The stoichiometry is fully flexible: transfers of carbon, nitrogen, phosphorus and silicate are tracked separately, with no fixed ratios, and the chlorophyll to carbon ratio varies depending on levels of light and nutrients. Kearney et al. (2021), reviewing the features of biogeochemical components in the next generation of global models, CMIP6, note that variable stoichiometry and inclusion of larger phytoplankton groups can improve the simulation of trophic transfers, while resolution of multiple detritus types enables fuller representation of nutrient circulation in complex, shallow-water environments. In this work, interactions at the seabed were represented by a simple remineralization scheme, where organic matter falling to the seabed is adsorbed and inorganic nutrients and carbon are returned to the water column at a rate proportional to their benthic concentration.

ERSEM was coupled online to POLCOMS, with the ecosystem model running in each model cell at each timestep. The model configuration and Southeast Asia domain are part of the Global Coastal Modelling System (Holt et al., 2009), which has previously been applied to investigate the effect of climate change on fisheries and marine food resources in many parts of the world (Blanchard et al., 2012; Barange et al., 2014; Fernandes et al., 2016).

The model domain covered the region from 99.1 to 139.0°E and 16.7°S to 23.9°N at a resolution of 0.1° latitude and longitude, with all sea areas within 200 km of the shelf break included ( Figure 1 ). Advective boundary conditions were applied at the open boundaries where the model domain is adjacent to the wider sea, enabling transfer of nutrients and other materials in and out of the domain. Model cells close to the open boundary have been removed from the analysis presented here to avoid boundary effects: these areas are shown as grey in Figure 1 .

Eight sub-regions were selected to sample the range of projected change across the model domain ( Figure 1 ). Four are coastal sites of key importance for biodiversity and sustainable development, as listed in the introduction. These are supplemented by four offshore areas: Box A is at 18-20°N, 114-116°E; Box B at 5.5-7.5°N, 109-111°E; Box C at 1.2°S-0.8°N, 125-127°E; Box D at 3.5-5.5°S, 127-129°E. These boxes were chosen to sample different parts of the region in places where in situ observations of nutrients and oxygen are available, and thus where the skill of the model in representing ocean processes could be assessed. Boxes A and B are both on the shelf edge, with waters 1000-3000m deep. Boxes C and D have variable bathymetry though without a strong shelf break; Box C has some shallower areas about 150 m deep but averages about 1000 m, most of Box D is more than 2000 m deep. All areas experience seasonally-reversing currents, weaker in Box C than for the other areas ( Figure 2 ).

Climate change was applied by using publicly available outputs from the Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor et al., 2011), a set of climate model experiments performed using protocols agreed by the World Climate Research Programme which inform the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC, 2021). Lateral boundary conditions were taken from the global model HadGEM2-ES (Jones et al., 2011, sourced from the Earth System Grid Federation archive, https://esgf-index1.ceda.ac.uk/projects/esgf-ceda/) and surface conditions from a regionally-downscaled atmospheric model driven by the same global model, HadGEM2-ES-RCA4, which was sourced from CORDEX Southeast Asia (http://www.ukm.my/seaclid-cordex/; Tangang et al., 2020). This pair of models was selected from the limited set available at the time the modelling was started, in 2019, because it provided both physical and biogeochemical boundary values in the ocean and a physically consistent regional-scale atmospheric component. The atmospheric model provided data at 0.22° resolution: temperature, pressure, wind and humidity at 6-hourly intervals and precipitation and radiation flux at daily frequency. The ocean model provided monthly values of temperature, salinity and current speed at 1° spatial resolution. It also provided annual concentrations of nitrate and inorganic carbon; these were then used to derive concentrations of phosphate and silicate and the monthly cycle of variation using ratios based on the present-day climatology from the World Ocean Atlas 2013 (Garcia et al., 2013a). The atmospheric and oceanic model values were used directly, without bias correction, to retain internal consistency. Atmospheric carbon dioxide levels for climate scenarios were set to the global levels defined for RCP4.5 and RCP8.5 (Meinshausen et al., 2011) for 2006 onwards and to the historical values used in CMIP5 for earlier years.

River inputs of fresh water and nutrients used average present-day values from the global model NEWS2 (Mayorga et al., 2010). Fixed values were used for nutrient concentrations; discharge values were adjusted to give a daily climatology based on historic flows reported by Dai et al. (2009). Future values for riverine discharge were based on this climatology, adjusted each year in line with domain-average precipitation. No information about future changes in river water quality was available, therefore concentrations of nutrients were kept constant at present-day values.

Initial conditions of temperature, salinity, oxygen and nutrients were taken from the World Ocean Atlas 2013 (Locarnini et al., 2013; Zweng et al., 2013; Garcia et al., 2013a; Garcia et al., 2013b) and total alkalinity and dissolved inorganic carbon from the GLODAPv2.2016b gridded dataset (Key et al., 2015; Lauvset et al., 2016). The model was run for an 11-year spin-up period using repeating yearly surface and boundary conditions which were the average for 1971-1980; at the end of the spin period conditions in both physical and ecosystem variables were stable. The model was then run for the period 1971-2098. For 1971-2005 the boundary and surface conditions came from the historical experiment of HadGEM2-ES and HadGEM2-ES-RCA4; for 2006-2098 the RCP4.5 and RCP8.5 experiments were used.

Model outputs were compared to observations to assess how well the model reproduces real-world conditions. Three-month mean conditions for sea surface temperature, currents and chlorophyll were plotted using values from the model and from satellite observations as listed in section 2.6; satellite values were masked to match the model domain, to make visual comparison easier. The period included was 1998-2017 in each case, to give a consistent 20-year period when satellite data were available. The modelled trends in sea surface temperature, surface chlorophyll, dissolved oxygen, net primary production and surface pH were compared to values from satellite observations and from the literature. The comparison period varied to use the longest observational dataset available in each case. The model trend was found by using a linear fit to the monthly mean values for the whole region. For all the comparison to observations, the model values for 2006-2017 were taken from the RCP8.5 run; the results using RCP4.5 were also calculated but were very similar and are not included here. Further comparison to in situ observations of temperature, salinity, nitrate, phosphate and oxygen from the World Ocean Database is provided in the supplementary material .

Present and future values for sea surface temperature and column total net primary production were compared using twenty-year averages for the present-day (2000-2019), mid-century (2040-2059) and end-century (2079-2098) and time series of annual mean values for 1980-2098. The present-day averages used RCP8.5 for 2006 onwards; tests with RCP4.5 gave very similar values. Annual mean values were calculated for the whole domain and for the sample areas shown in Figure 1 ; anomalies were calculated from these as the difference between the annual value and the mean for 1980-2005, i.e. the historical period as defined for CMIP5.

Mid-century and end-century change in a wider range of variables was summarized by calculating:

The index of change enables the change in different variables to be compared to each other and to the range of conditions experienced in the present day: change that goes significantly outside the current range is likely to put more stress on organisms and ecosystems. The variables included were selected to give a range of ecosystem-relevant parameters: surface and bottom temperature, surface salinity, column total net primary production, phytoplankton biomass and zooplankton biomass and surface and bottom level oxygen and pH.

It was only possible to downscale one model in this study, but in order to put the outputs in the context of other modelling work they have been compared to a range of global model outputs from the CMIP5 set, including the model that is the parent for the regional modelling, HadGEM2-ES. The global model outputs were interpolated to the model grid using a nearest-neighbor method, then masks were applied to select the sample areas shown in Figure 1 . Mean values for each area were calculated and compared to regional values and to satellite; these are shown as anomalies in section 3 and as absolute values in the Supplementary Material . The Supplementary Material also includes some comparison of the global model outputs to satellite and in situ observations.

Satellite-derived products sourced from the Copernicus Marine Service (https://marine.copernicus.eu) were used to make plots comparing the model outputs to observations:

The same sea surface temperature and chlorophyll products were used for quantitative assessment of trends over time ( Table 1 ). Information about trends in dissolved oxygen, net primary production and surface pH was derived from the literature and sources are given in Table 1 .

CMIP5 global model outputs were sourced from the Earth System Grid Federation portal (https://esgf-index1.ceda.ac.uk/projects/esgf-ceda/).

Three-month average modelled surface temperature, currents and chlorophyll for 1998-2017 compared to satellite-derived values are shown in Figure 2 . The mean and seasonal variation in sea surface temperature is generally captured well by the model, the exceptions are that April-June (AMJ) temperatures are higher than observed in the west and July-September (JAS) temperatures are lower than observed in the north ( Figure 2A ). Surface current speeds are higher than satellite derived (geostrophic + Ekman) currents, but the observed circulation patterns are captured by the model, including seasonal reversal due to the monsoon climate in the north of the region ( Figure 2B ). The range of surface chlorophyll is much greater than observed by satellite, with higher highs and lower lows ( Figure 2C ). However, the spatial and temporal variations are similar, with the highest values along coasts, especially in the southwest and southeast, and the higher chlorophyll months being October to March in the northern part of the region and July to September in the south.

Compared to the parent global model, chlorophyll estimates by the regional model tend to be closer to satellite values at the coastal sites but not in open waters ( Figure S2 ).

Table 1 shows long-term trends from the model compared to values from satellite observations and reported in the literature. The model outputs and satellite observations both show a rising trend in sea surface temperature for the period 1985-2019, though the observed trend is smaller than modelled (see also Figure 3 ). The magnitude of change in surface chlorophyll concentration for 1998-2017 is small in both model outputs and satellite observations. Modelled trends in dissolved oxygen concentration are smaller than observed but have the same spatial distribution, with increases in the north and west and decreases in the south and east. The model outputs show little or no change in net primary production since 1998, which is consistent with the findings for this region in a global analysis of satellite and in situ observations (Kulk et al., 2020). The model shows a small decrease in surface pH, -0.01 to -0.02 pH units per decade, except for shallow areas in the northwest of the region where there is a small increase. There are few measurements of pH for the Southeast Asia region to compare this with, but the IPCC reports a global decrease of 0.013 to 0.03 pH units per decade for 1990-2014 (IPCC, 2019).

Further comparison to in situ observations is given in the supplementary material . Overall, the model outputs show a reasonable match to the range and median of observations ( Figure S3 ). Agreement is worst for Box A, where temperatures and oxygen are underestimated and nitrate is overestimated. Nutrients are also overestimated in Box C.

This section presents model outputs for conditions to the end of the 21st century under the moderate RCP4.5 scenario (atmospheric greenhouse gas concentrations rising until mid-century and then stabilizing) and the more extreme RCP8.5 (greenhouse gases continuing to rise throughout the century). Sea surface temperature and column total net primary production are discussed in some detail; change in other environmental variables is summarized and discussed more briefly.

Warming seas mean that by the end of this century some parts of the Southeast Asian seas may experience average temperatures not seen in the region at present ( Figure 4 ). Species that live at or near the sea floor are particularly likely to experience conditions of temperature, oxygen and pH outside the current range of variability, to which they may not be adapted ( Figure 9 ). As a response to rising temperatures, some species populations may be able to move with the present-day temperature contours, seeking to maintain optimum conditions for their growth, reproduction and survival (Pinsky et al., 2013; Poloczanska et al., 2013). However, this adaptation strategy has limited success near warmer regions of species distributions, and does not apply to all types of species, especially sessile organisms – for example, range expansion of corals is restricted because they rely on the tolerance range of their algal symbionts to adapt to new environments (Tonk et al., 2013), and these symbiotic algae are tightly dependent on light and temperature near the sea surface.

All the environmental changes mentioned throughout this article will have consequences for the fisheries sector. This is critical in a region where there are already significant challenges due to overcapacity and declining catches, exacerbated by poorly regulated fisheries (Pomeroy, 2012; Teh et al., 2017). Further modelling is needed to assess the implications of climate change for different species, and the projections presented here could provide input data for such modelling. We can expect that populations of reef-dwelling species such as grouper are likely to decline and small pelagics such as sardine, clupeides and mackerel will re-distribute in response to changes in habitat conditions. Therefore fishers may need to travel further to find their target species, shift to catching different species, perhaps requiring an investment in alternative gear or larger vessels, or in the worst cases be faced with declining catches and no incoming new species (Barange et al., 2018). This is especially a concern in equatorial regions, such as Southeast Asia, where declining local diversity resulting from poleward retreat of distribution leading edges is not necessarily compensated by new species arrivals. Therefore, adapting to climate change will be a significant challenge, especially for the many small-scale fishers in the region, and good management will be essential for protecting livelihoods (Lam et al., 2020). There is a need for modelling of future fish productivity and distribution to inform decision-making.

Increased temperature and reduced salinity, as seen in these projections, may result in increased incidence of harmful algal blooms (HABs). The consequences of such HAB events include reduced water quality and toxin build-up in fish and shellfish, with potential subsequent impacts on human health and food security (GEOHAB, 2010; Gobler, 2020; Young et al., 2020; Trottet et al., 2022). Rising temperatures can also increase the risk of disease in cultured fish (Reverter et al., 2020), shellfish (Allison et al., 2011) and seaweed (Largo et al., 2017), and can lead to a rising number of jellyfish blooms or invasions that may affect aquaculture (Xu et al., 2013; Bosch-Belmar et al., 2020). Changes in the timing of seasons are also likely to affect aquaculture production cycles, where activities are tightly timed to sharp climatic differences between monsoon/inter-monsoon periods. Such changes can be seen in the projected changes for coastal areas, notably Palawan and Taka Bonerate-Kepulauan Selayar ( Figure 5 ). A reduction in the predictability of seasonal cycles often leads to reduced harvests (Handisyde et al., 2006; Hamdan et al., 2015). The projections presented here could be a useful tool for assessing the impacts of climate change on phenology and future potential for aquaculture in different parts of the region.

All the analyzed coastal areas have projected end-century temperature increases close to 1.5°C under RCP4.5, enough to cause thermal stress leading to coral bleaching, while the 2.7°C increase seen under RCP8.5 is likely to cause widespread loss of coral (Lough et al., 2018). The smallest temperature increases were seen at Cu Lao Cham, but at 1.5°C (RCP4.5) or 2.5°C (RCP8.5) even these are too high to prevent coral damage (Dao et al., 2021). Ocean acidification and the overall alteration of the ocean carbonate system resulting from rising atmospheric CO₂ levels creates further stress to coral reefs by affecting the ability of reef organisms to maintain sufficient calcification rates in the face of increased dissolution rates and, in extreme cases, prevent the deposition of carbonate minerals needed for skeleton construction through the occurrence of insufficient saturation levels (Eyre et al., 2018). Our projections show surface pH decreasing, though not beyond the range currently experienced; changes in bottom level pH go outside the current range at many locations. Any acidification will act as an additional stressor on coral reefs and make it more difficult for them to recover from bleaching events. The conditions under which coral reefs are able to recover can be complex (Graham et al., 2015) and this has not been considered in the current study.

We were only able to downscale one global model in this work, because of resource constraints and the limited availability of downscaled atmospheric projections at the time the work was started. More regional model projections, based on a range of recent global climate models, are needed to give information about uncertainty. These should follow best practice methods developed from downscaling other regions (Drenkard et al., 2021) including a careful choice of global models to span a range of possible futures, attention to bias correction, and using time periods long enough to capture long-term trends. A downscaled atmospheric model should be used at the surface to capture dynamic features at scales relevant to ecosystem processes; for this region that includes the intense storms associated with tropical cyclones, which are captured better by regional than global models, and monsoon winds influencing coastal upwelling. More regional atmospheric models are now available (Tangang et al., 2020) so it is possible to select from multiple options which have both downscaled atmospheric projections and ocean biogeochemistry from the same global model; in due course this will include downscaled versions of the most recent generation of global models, CMIP6.

A range of regionally-appropriate biogeochemical models should also be used within the regional modelling system. The global model outputs we looked at vary widely in their biogeochemical outputs ( Figure 7J , Figure S5 ), indicating that it would be valuable to include a number of models to get an understanding of the inter-model uncertainty. CMIP6 models include biogeochemical components with greater complexity (Kearney et al., 2021) which could potentially represent the region better, and it would be useful to explore their range of response for Southeast Asia and compare them to CMIP5. But it will also be important to explore the outputs of different models used within a regional modeling system that can simulate the shallow shelf waters and steep slopes of the seas in this region.

Future regional modelling work can build on the results reported here. A useful next step would be to investigate the reasons for the over-prediction of nutrients in the north of the region. Possible causes include inaccurate initial conditions, giving persistently high subsurface nutrients available to be mixed into the eutrophic layer ( Figure S3 ), an over-supply of river-borne nutrients, or overly strong surface currents ( Figure 2 ) leading to too much mixing of nutrients from below. It would also be useful to have information about projected future change in river water quality, especially for large rivers such as the Mekong. This is likely to have high uncertainty, but the sensitivity of the system to river inputs can still be explored.

A further aspect of regional modelling that would benefit from improvement is future changes in storminess, which is of high concern for people in the region. Typhoons disrupt ocean ecosystems through water column mixing by strong winds and the freshening effect of heavy rain. They can also cause flooding and surging which may damage fishing gear and cages and cause long-term damage to coral reefs and other habitats (Harmelin-Vivien, 1994; Latypov and Selin, 2012; Safuan et al., 2020). Increased typhoon activity leads to beach erosion, causing damage to infrastructure and property for communities living close to the shore (Shariful et al., 2020). All of these effects are exacerbated by sea level rise, which poses an additional threat to the large coastal population of this region (Rowley et al., 2007; Neumann et al., 2015; Saleh and Jolis, 2018). There are indications that climate change is resulting in the increased intensity and frequency of typhoons and flooding in SE Asia (Loo et al., 2015), although considerable uncertainty remains about both historical and future trends (Lee et al., 2012; Ying et al., 2012; Knutson et al., 2020). Regional models can simulate stronger storms than global models, because they have higher resolution, but the uncertainty in the projections remains high. We investigated the surface wind speeds in the regional atmospheric model HadGEM2-ES-RCA4, which was used as input to the marine model reported here, for any changes in the strength, frequency or timing of storms. There was some indication of an increase in the number of days per year with strong winds, but not in maximum wind strength or in the pattern across the year. By contrast, Herrmann et al. (2020), based on a much more thorough analysis of outputs from a similar regional model (CNRM-CM5_RegCM4), found a decrease in projected wind speeds, in most of the region and all seasons, except for some increase in average speeds for December to February in the north of the region. The number of tropical cyclones also decreased in all seasons. More work is needed to resolve this uncertainty.