The Spatial and Temporal Variability of Air-Sea CO2 Fluxes and the Effect of Net Coral Reef Calcification in the Indonesian Seas: A Numerical Sensitivity Study

A numerical model system was developed and applied to simulate air-sea fluxes of CO2 and coral reef calcification in the Indonesian Seas and adjacent ocean basin for the period 1960 to 2014 on a fine resolution grid (ca. 11 km) in order to study their response to rising sea water temperatures and CO2 concentrations in the atmosphere. Results were analyzed for different sub-regions on the Sunda Shelf (Gulf of Thailand, Malacca Strait, Java Sea) and show realistic and different levels, signs and pronounced temporal variability in air-sea CO2 flux. The Gulf of Thailand changes from an atmospheric CO2 sink during the boreal winter to a CO2 source in summer due to higher water temperatures, while other sub-regions as well as the entire averaged Sunda Shelf act as a continuous source of CO2 for the atmosphere. However, increasing atmospheric CO2 concentrations weakened this source function during the simulation period. In 2007, the model simulations showed even a first flux inversion, in course of which the Java Sea took up CO2 . The simulated trends suggest that the entire Sunda Shelf will turn into a permanent sink for atmospheric CO2 within the next 30 to 35 years if current trends remain constant. Considering the period between 2010 and 2014, coral reef calcification enhanced the average CO2 emission of the Sunda Shelf by more than 10 % from 15 to 17 Tg C yr−1 due to lowering the pH and increasing the partial pressure of CO2 in surface water. During the entire period of simulation, net reef calcification decreased although increasing seawater temperature mitigated effects of reduced CO2 emission and the resulting decrease of the pH values on reef calcification. Our realistic simulation results already without consideration of any biological processes suggest that biological processes taking up and releasing CO 2 are currently well balanced in these tropical regions. However, the counteracting effects of climate change on the reef calcification, on other biological processes and the carbonate system need to be investigated in more detail. SST ◦increased by about 0.6 C during the last 55 years, while SSS decreased by about 0.7 psu.

A numerical model system was developed and applied to simulate air-sea fluxes of CO 2 and coral reef calcification in the Indonesian Seas and adjacent ocean basin for the period 1960-2014 on a fine resolution grid (ca. 11 km) in order to study their response to rising sea water temperatures and CO 2 concentrations in the atmosphere. Results were analyzed for different sub-regions on the Sunda Shelf (Gulf of Thailand, Malacca Strait, Java Sea) and show realistic and different levels, signs and pronounced temporal variability in air-sea CO 2 flux. The Gulf of Thailand changes from an atmospheric CO 2 sink during the boreal winter to a CO 2 source in summer due to higher water temperatures, while other sub-regions as well as the entire averaged Sunda Shelf act as a continuous source of CO 2 for the atmosphere. However, increasing atmospheric CO 2 concentrations weakened this source function during the simulation period. In 2007, the model simulations showed even a first flux inversion, in course of which the Java Sea took up CO 2 . The simulated trends suggest that the entire Sunda Shelf will turn into a permanent sink for atmospheric CO 2 within the next 30-35 years if current trends remain constant. Considering the period between 2010 and 2014, coral reef calcification enhanced the average CO 2 emission of the Sunda Shelf by more than 10% from 15 to 17 Tg C yr −1 due to lowering the pH and increasing the partial pressure of CO 2 in surface water. During the entire period of simulation, net reef calcification decreased although increasing seawater temperature mitigated effects of reduced CO 2 emission and the resulting decrease of the pH values on reef calcification. Our realistic simulation results already without consideration of any biological processes suggest that biological processes taking up and releasing CO 2 are currently well balanced in these tropical regions. However, the counteracting effects of climate change on the reef calcification, on other biological processes and the carbonate system need to be investigated in more detail. SST increased by about 0.6 • C during the last 55 years, while SSS decreased by about 0.7 psu.

INTRODUCTION
The ocean-atmosphere CO 2 exchange process is quite well known (Wanninkhof, 1992, 2014, andothers). It depends on the concentrations or partial pressures of CO 2 in the ocean surface water and the atmosphere and the gas transfer velocity. The oceanic CO 2 concentration is influenced by its solubility and, therefore, by the water temperature and salinity as well as the marine buffer system (alkalinity). It can change temporally and spatially due to physical and biogeochemical processes. The transfer velocity is related mainly to the wind speed above the ocean surface and the Schmidt number and temperature of the ocean surface water. In general, an atmospheric increase of CO 2 content leads to less CO 2 transfer from the ocean to the atmosphere or to more transfer into the other direction (or to a flux direction change). With this, the CO 2 concentration in the ocean increases resulting in the well known ocean acidification with all consequences to the marine ecosystem.
Air-sea fluxes of CO 2 between continental shelf and the atmosphere are still poorly quantified for the Indonesian Seas due to the poor data coverage, the complexity of the interacting processes and their response to increasing CO 2 concentrations in the atmosphere and the associated global warming (Cai et al., 2006;Chen and Borges, 2009;Doney et al., 2009;Takahashi et al., 2009;Borges, 2011;Wang et al., 2012;Bauer et al., 2013;Jones et al., 2014). Especially in the Indonesian Seas, independent model based estimated CO 2 fluxes along the coast disagree with air-sea fluxes of CO 2 derived from global ocean circulation models (Laruelle et al., 2014). The Sunda Shelf covers more than 20% of the global tropical shelf seas defined as regions less than 200 m deep between 23.4 • southern and northern latitudes, i.e., 1.89 out of 9.21 10 6 km 2 , according to the 30 ′′ × 30 ′′ digital topography we used for this work (Farr et al., 2007). The Indonesian Seas are furthermore hot spots of coral diversity suffering from over-exploitation, global warming, and ocean acidification like all corals world wide (Cinner et al., 2016). The total area of corals in this region, 127,028 km 2 , covers almost 30% of the global warm water coral reefs based on the global map from UNEP-WCMC et al. (2010) (423,582 km 2 ). To our knowledge, Kartadikaria et al. (2015) are the first and so far only authors, who compiled detailed data on air-sea fluxes of CO 2 from the Indonesian Seas.
In contrast, the complexity of coral reef calcification, the interaction with reef pore water fluxes and changes in seawater chemistry have been studied intensively (Lough and Barnes, 1997;Langdon et al., 2000;Kleypas et al., 2006;Hoegh-Guldberg et al., 2007;Silverman et al., 2007;McGowan et al., 2016;Page et al., 2017, and many others). The application of numerical models incorporating observational findings and studying the impacts of rising CO 2 concentrations and sea water temperatures on the air-sea exchange of CO 2 and the reef calcification is missing so far on a horizontally fine resolved grid. Therefore, we developed a warm water coral reef calcification module for the regional biogeochemical ecosystem model ECOHAM and applied this on a fine-resolution 3D grid to the region of the entire Indonesian Seas to simulate the chemical carbonate system in this tropical marine environment.
The overarching goal is to get a principal idea about the spatial distribution and the long-term development of geochemical parameters, especially the exchange of CO 2 between ocean and atmosphere, under changing physical conditions (climate change and increasing atmospheric CO 2 content) and under the influence of carbonate precipitation with reliable rates from realistically distributed coral reefs. For this, a numerical model experiment has been performed to simulate the carbonate system in the tropical ocean (1) entirely without biological processes but under the influence of climate change and (2) additionally with aragonite precipitation only from coral reefs. The exclusion of biological processes in part (1) of our model experiment enables us to investigate the dependency of the carbonate system solely on the physical environment and on the changing atmospheric CO 2 content. The consideration of a net coral reef calcification rate in part (2) of our experiment helps to investigate the undisturbed effect of net coral reef calcification on that carbonate system and vice versa the effect of a changing physical environment on the calcification rate.
In the next section, we introduce our setting of the numerical models applied to simulate the hydrodynamics and to perform the experiment for the carbonate system without and with calcification from coral reefs. Then, hydrodynamic and the ecosystem model results will be verified. Thereafter, ecosystem model results will be presented as long-term time series for different sub-regions as well as for certain locations with different coral reef coverage. A discussion and a summary follow the data presentation.

MATERIAL AND METHODS: APPLIED MODEL SYSTEM
The applied numerical model system consists in principle of a global general circulation model, which provides forcing data (sea surface height, temperature, salinity) at the open boundaries of a regional circulation model. Furthermore, a global hydrological model delivered river discharge data to the regional circulation model. This setup aimed at obtaining most realistic 2D and 3D time series for the physical conditions of the Indonesian Seas, especially for the water mass transports (velocities) and the 3D temperature and salinity distributions. Subsequently, an ecosystem model was applied on the same grid as the regional circulation model, of which the aforementioned physical parameters were used to simulate the chemical carbonate system.
Our regional model domain is presented in Figure 1. The horizontal grid resolution amounted to 6 min (approximately 11 km). In the vertical direction, 39 z layers with increasing thicknesses from 5 m at the surface to 550 m at greater depths resolved the bathymetry, which was limited to a maximum depth of 6,260 m.

Hydrodynamic Model
The global circulation model MPI-OM (Max Planck Institute Ocean Model, Marsland et al., 2003;Jungclaus et al., 2006), which is the oceanic part of the climate model system used for the IPCC world's climate simulations, was applied to the global ocean FIGURE 1 | Regional hydrodynamical model domain with bathymetry. Red dots: river input points; X: locations to show single point time series of simulated carbonate system data. Isobath of 200 m is also drawn to depict the Sunda Shelf limit. Along the upper and right axes, model grid indices are shown. on the TP04 grid (Jungclaus et al., 2013), having a horizontal resolution of 24 min (approximately 44 km at lower latitudes) and 40 z layers of increasing thicknesses from top to bottom. Meteorological forcing data were taken from NCEP/NCAR (Kalnay et al., 1996). The reason for this first application step was to obtain high quality forcing data for the subsequently applied regional ocean circulation model HAMSOM (Hamburg Shelf Ocean Model, Backhaus, 1985;Pohlmann, 1996Pohlmann, , 2006Mayer et al., 2010), which used the vertical sections of temperature and salinity and the sea surface height as six hourly time series for its lateral open boundaries to the Indian Ocean, Pacific Ocean and South China Sea (Figure 1). HAMSOM is a prognostic baroclinic regional circulation model solving the momentum and tracer transport equations semi-implicitly on the Arakawa-C grid. Meteorological forcing data were taken from the same source (NCEP/NCAR) as for the global model.
Additionally, HAMSOM considers freshwater input through rivers into the ocean. These river discharge data were generated with the MPI-HM (Max Planck Institute Hydrology Model, Stacke and Hagemann, 2012), which is a global hydrological model solving the water balance on the land surface. The MPI-HM regularly takes part in model inter-comparison projects (Haddeland et al., 2011;Schewe et al., 2014). The input points for freshwater in the regional model domain are depicted by red dots in Figure 1.
For this investigation, long-term simulations were performed to analyse the development of properties over a period of more than 60 years. Therefore, daily averaged values were used to force the subsequently applied ecosystem model simulating the carbonate system.

Ecosystem Model
For the simulation of the carbonate system in the Indonesian Seas, ECOHAM (ecological model, Hamburg, Pätsch and Kühn, 2008;Kühn et al., 2010;Lorkowski et al., 2012;Große et al., 2016) has been applied on the same grid as the circulation model HAMSOM. ECOHAM is an NPZD-type ecosystem model simulating the cycles of nutrients (N, P, C, Si), of all species of dissolved inorganic carbonate (DIC), of O 2 and of total alkalinity (TA) and was originally developed to simulate biogeochemical processes including phyto-and zooplankton dynamics in the North Sea, which is located in the temperate climate zone. Parameters and limits in the process formulations of this model, in particular for biological processes, are adapted to that climate zone but not to the tropical zone. As a first step toward a fine resolution simulation of the carbonate system in the tropical ocean, we switched off the biological components in ECOHAM and used only the carbonate system part of the model. With this, it is possible to investigate the sensitivity of the carbonate system purely on physical environmental parameters plus the CO 2 concentration in the air. In a second step, we added a new module to simulate the CaCO 3 (only aragonite) precipitation due to calcifying corals to investigate also their influence on the carbonate system. This enables us to investigate the undisturbed effect of coral reef calcification on the carbonate system. Nutrients and their dynamics were not considered in these simulations (no biological components).
The principle way of simulating the carbonate system by ECOHAM is the following: transport and redistribution of concentrations of TA and DIC through the model grid by advection and diffusion with the velocities and turbulent diffusion coefficients provided by the regional circulation model HAMSOM. From the concentrations of DIC, the inorganic carbonate species dissolved CO 2 , carbonic acid, bicarbonate and carbonate ions are calculated according to their solubility constants. Due to the CO 2 flux between ocean and atmosphere, DIC concentration can increase or decrease in the upper model layer. To provide the anthropogenic effect, a formula calculated the time dependent atmospheric CO 2 concentrations with best fit to observational data (Figure 2A), which have been kindly made available by Dlugokencky et al. (2017).
As initial state and for the lateral open boundaries, different data sets of observed TA and DIC concentrations were taken from the CDIAC database Carbon Dioxide Information Analysis Center, RRID:SCR_006999, Johnson et al., 2002;Feely et al., 2007;Hydes et al., 2012;Kamiya et al., 2012;Tilbrook and Wijffels, 2012). These data from different ship cruises are very sparse in space and time for our region of interest. They were spatially extra-and interpolated to obtain a 3D distribution of DIC and TA concentrations in the model domain, which is not synoptic. However, there is currently no other way of creating such a fine resolved distribution of these parameters for this large domain.

Module to Simulate Net Coral Reef Calcification
For the simulation of precipitation of CaCO 3 (aragonite) from coral reef calcification, we basically follow the equation for the rate of calcification (Burton and Walter, 1987) and assume that the saturation of aragonite is drastically increased in the calcifying fluid due to the metabolism of the corals according to McCulloch et al. (2012).
G is the calcification rate (µmol m −2 h −1 ), a the saturation state of aragonite, k the rate constant (µmol m −2 h −1 ), n the order of reaction (Burton and Walter, 1987) and k a the solubility constant for aragonite. k and n are functions of water temperature according to McCulloch et al. (2012), k a is a function of temperature, salinity and pressure (Mucci, 1983;Millero, 1995). Calcification is supposed to produce additional CO 2 in the ocean water: Biomineralization is still a poorly understood process. Organic matter influences crystal growth and the morphology of carbonate skeleton, but calcification effects alkalinity and DIC as described in Equation. 3 according to our current understanding, as pointed out by Tambutté  The calculated amounts of carbonate ions, which are taken out of the water due to the "carbonate precipitation, " are then deducted correspondingly from the DIC and TA pools.
Coral growth reduction and/or mortality during extreme events as seen during the ENSO periods in 1997/8 and 2015/16 are not considered in our simulations. Accordingly, increasing temperatures enhance the calcification process during the ENSO events in 1997/1998, when thermal stress caused an abnormally high coral mortality in the Indian Ocean (Wilkinson et al., 1999). Simulated coral growth rates are thus considered with great care. Due the short duration of extreme events, their overall impact on reef calcification and air-sea fluxes of CO 2 are assumed to be negligible, even though large areas of mortified coral reefs do have an impact on the local marine system.
In our model simulations, all the above described calcification related processes take place only in grid cells with a certain reef coverage. Using the UNEP global map of warm water coral reefs (UNEP-WCMC et al., 2010), which gives the presence of coral reefs on a grid with 30 × 30 s mesh size, a distribution of corals in our domain with 6 × 6 min resolution was calculated as coverage in percent of the horizontal area of each water column (see Figure 2B). The vertical location of the corals in the 3D model grid is the bottom layer but no deeper than 22 m (lower end of third layer).

Validation of Model Results
For the investigation of carbonate system related processes, it is essential to use a reliable distribution of temperature and salinity as well as transport of water masses. As examples to prove the high quality of our model results, the following sections will show the agreements of simulated and observation derived transports through key sections in our model domain to prove that the general simulated circulation is close to reality. Furthermore, a comparison between sea surface temperature and salinity data at certain locations will also prove realistic physical model results.

Simulated Averaged Transports Through Selected Sections
The averaged simulated transports through some important sections in the Indonesian Seas for the periods 2004-2006 and 1970-2006 are presented in  From the transport distributions, we can detect a principle agreement of the simulated with the expected circulation in the Indonesian Seas according to the cited literature. Remarkable differences occur for the Ombai Strait and Timor Passage transports, which is probably due to the fact that the southern part of the eastern open boundary cuts directly through the Arafura Sea between New Guinea and Australia. Obviously, in the model, a part of the water entering from the Pacific Ocean leaves through this open boundary instead of recirculating in the Arafura Sea and subsequently leaving the area through Ombai Strait and Timor Passage into the Indian Ocean.
The Lifamatola Strait value of −2.5 Sv (van Aken et al., 2009) represents only the depth range below 1250 m. The authors guess that there is rather a northward transport in upper regions of the passage leading to a total southward transport through this section of even less than 2.5 Sv. Our results of −3.8 Sv, however, correspond to the second cited value, −3.2 Sv, from Feng et al. (2013).

Comparison of Simulated SST and SSS With
Satellite Data The regional model HAMSOM agrees also well with real conditions concerning sea surface salinity (SSS) and temperature MODIS (Moderate-resolution Imaging Spectroradiometer) data from both EOS-Aqua and EOS-Terra satellites have been averaged temporally and horizontally over the corresponding depicted regions to obtain regional time series with monthly values. MODIS data have an accuracy of 0.4 • C and represent surface squares of 4.6 km × 4.6 km. Simulated SST time series were created by averaging the upper layer temperature for the same regions and then temporally from daily to monthly data. When comparing the time series, we have to keep in mind that MODIS data are "skin temperatures" (Brown and Minnett, 1999). They are influenced by diurnal solar surface heating and by slight cooling due to evaporation (Minnett, 2003) with a heating effect much stronger than the cooling effect. An average of day and night remote sensing data leads therefore to slightly higher temperature data than the surface water would show. In contrast, the model results are "bulk temperatures" for the model's surface layer with a thickness of about 5 m plus the model's water level.
The comparison displayed in Figure 4 clearly shows that the seasonality as well as the level are in very good agreement, the seasonal peak values are over-or underestimated by no more than 1 K. It has been explained in the previous paragraph that the resulting MODIS averages are slightly higher than subsurface temperature. Therefore, model data are expected to be lower than these remote sensing data.
For SSS comparison, SMOS (soil moisture and ocean salinity sensor) derived data represent-with an accuracy of 0.4 psuthe upper ≈ 1 cm of the ocean surface (Boutin et al., 2014). They are available for horizontal areas of 1 • × 1 • , averaged from 3 days maps to monthly means. The simulated SSS data are again the values of the upper model layer (thickness: 5 m plus water level). The data were taken at the positions indicated in the map of Figure 4 and averaged to monthly values. The comparing time series are presented in panels 1-7 of Figure 4. Also here, a very good agreement of simulated and observed data is obvious for both the level and the seasonal variation if there is any.
Differences are usually not more than 1 psu. It is impressive to see the seasonal variation at location 7, eastern Java Sea: During the first and last quarters of a year, the location is effected by relatively fresher water masses coming from the South China Sea, whereas during boreal summer season, this region is effected by the saltier Makassar Current coming from the North Pacific Ocean (Mayer and Damm, 2012).
A slight overestimation of SSS from the model is expected because of the thickness of the upper layer. Precipitation has a much larger impact on the remote sensing data than on the model data. Since the model does not resolve any vertical profile within the upper 5 m of the ocean, fresh water input from precipitation is immediately mixed into the upper 5 m layer, while the satellite sensor detects the fresher water within the upper ≈ 1 cm of the ocean. Because our locations are in the tropical zone with high rain fall rates, we expect the remote sensing data to be rather lower than the model data, which can indeed be seen in the comparing time series of Figure 4.

Comparison of Simulated and Observed Carbonate System Data and Calcification Rates
During a research cruise (MTK-2012, October 2012) with the vessel MatahariKu, the surface water pH and partial pressure of CO 2 in the water (pCO 2 ) was continuously measured by an underway pCO 2 system (SUNDANS) and a ferry box in the Bantan Bay (5.966 • S, 106.134 • E), a near-shore site about 70 km west of Jakarta. Wit et al. (2015) describe the method in detail, and Figure 5 presents the data. Based on the pH and pCO 2 values, the CO2sys program was used to calculate DIC and TA. These four parameters were compared with simulated data. The simulated data present daily means for each of the upper three layers at the time of the cruise for the grid point at the location of the measurements in the Banten Bay.
This comparison shows that the model simulation overestimates TA by about 100 µmol kg −1 and pH value by ca. 0.15 and underestimates pCO 2 . Considering that Banten Bay is a near-shore and polluted area, the deviation might reflect the fact that the model TA is rather an open ocean value and not effected by near-shore ecosystems, since their benthic and pelagic biochemical activities are not considered in the simulations-e.g., Borges et al. (2005) show that the carbonate system in the coastal ocean exhibits different behavior compared to the open ocean (for example different air-sea CO 2 flux). Another possible reason for the differences is the quality of input data for the model, which is the best we could achieve within our study: TA is an averaged conservative tracer, interpolated from the CDIAC database and prescribed at the open model boundaries and transported by the currents into the model domain. Thus, it might not capture the specific situation in the Banten Bay.
A comparison of simulated carbonate system parameters (values and trends) with published estimations will be presented in the next section "Results", also the good agreement of simulated calcification rates with published estimations from observations. Figure 6 presents the horizontal distribution of CO 2 flux between ocean and atmosphere in two different seasons (19 January and 13 July 2014). The overall distribution of uptake and outgassing is in good agreement with Kartadikaria et al. (2015), who analyzed a collection of observational data from the last few decades. Also, our simulated increase in pCO 2 of 76-120 µatm century −1 is within their range of 60-380 µatm century −1 . In contrast to the Gulf of Thailand, parts of the Malacca Strait and the Karimata Strait, where the air-sea fluxes of CO 2 reverse on a seasonal cycle, the entire Sunda Shelf acts as a continuous source of CO 2 to the atmosphere. Consideration of coral reef calcification in the model (right panels) hardly affects the general pattern but causes small scale local differences in reef dominated regions. In these areas (e.g., northwest of Bangka Strait, in the Karimata Strait next to Pulau Belitung, and in the Spermonde Archipelago), calcification can even change the sign of the air-sea fluxes and convert CO 2 sinks into CO 2 sources for the atmosphere.

RESULTS
To obtain spatial averages for the different sub-domains, 3D variables (TA, DIC, pCO 2 , and pH) were volume weighted averaged, 2D variables (atm-oce CO 2 flux, SSS, and SST) were area weighted averaged, and directly reef dependent variables ( a and calcification rate) were area weighted averaged over the corresponding coral reef areas according to the map in Figure 2B.

Carbonate System in Different Regions
The following carbonate system relevant parameters are presented in this section for the period 1960-2014 as spatial averages for different sub-domains: TA, DIC, pCO 2 , air-sea CO 2 flux, pH value in seawater scale, a in seawater and in the calcifying fluid, aragonite calcification rate, SSS and SST. Figure 7 shows time series of yearly running means of the simulation results with corals for the inner Sunda Shelf, Gulf of Thailand, Malacca Strait and Java Sea, while Figure 8 displays time series of monthly values from both simulations without and with corals for the Java Sea (as an example) to display the role of coral reef calcification and the seasonality. In Table 2, average values for the last 5 years of simulation with coral reef calcification and trends in the time series from the last 50 years of simulation are given for all parameters and all sub-domains.
Looking at the long-term time series of the parameters in Figure 7, we see increasing values for DIC, pCO 2 , and SST, and decreasing values of pH and a both in seawater and in the calcifying fluid, while the calcification rates and the CO 2 emissions slightly decrease.
The simulated range of pCO 2 increase of 76-120 µatm per century, corresponds well with values listed by Takahashi et al. (2009,  The simulated CO 2 flux to the atmosphere in the Sunda Shelf is in good agreement with estimations of Bauer et al. (2013). Our contemporary sea to air flux comes to ∼8.4 g C m −2 y −1 on the average and about 0.017 Pg C y −1 in total for the inner Sunda Shelf area ( Table 2). This is ∼21% of the global low latitude continental shelves' contribution as published by Bauer et al. (2013, 0.08 Pg C y −1 ). Our inner Sunda Shelf covers about 23% of the global low latitude (< 30 • ) continental shelf area (1.89 10 6 km 2 out of 7.91 10 6 km 2 FIGURE 5 | Simulated and measured TA (alk), DIC, partial CO 2 pressure (pco2), and pH in the surface water at Banten Bay, just north of western tip of Java. For the measured data, TA and DIC were calculated from the other two parameters. Labels along the x axis are Julian days of 2012 (274: 30 Sep).
according to Bauer et al. (2013). Laruelle et al. (2014) estimate a sea to air flux of 0.023 Pg C y −1 for their region "South East Asia" representing an area of 2.3 10 6 km 2 , which is similar to our results. Cai et al. (2006) mention a slightly higher outgassing rate of 12 ± 5 g C m −2 y −1 for low latitude western and eastern boundary current shelves. Assuming that their estimations are based on older data and that there is a trend in the flux toward decreasing CO 2 emissions, our model results match their findings. Chen and Borges (2009) report in their Table 3 a much lower CO 2 outgassing rate of 1.2 g C m −2 y −1 for tropical continental shelves in the eastern Indian and western Pacific Oceans.
From the SST time series, we can observe temperature peaks in nearly all sub-domains for the year following a strong El Niño event (first half of 1973, 1983, 1998, 2010), the effect of which can be clearly detected in the calcification rate showing peaks at the same times, most obvious for the Gulf of Thailand, where the pH value is highest. The SSS time series do not show any ENSO response-except for the Java Sea, where we get a salinity peak during strong El Niño phases reflecting reduced precipitation rates.
To get an idea about the seasonality and the effect of calcification on the above mentioned parameters, Figure 8 presents the same parameters as monthly averages only for the Java Sea but as result from the simulation without biology and with calcification. There are clearly two different seasonal signals visible in the time series, which are introduced by the altering physical environment: Due to the Java Sea's geographical location close to the equator, the sun passes the region's zenith twice a year, producing two SST peaks per year (usually in Apr and Nov), with consequences on the solubility of CO 2 in the seawater and subsequently on pCO 2 and other parameters. Additionally, the predominant water flow direction alters seasonally: during boreal summer monsoon season, the westward flow transports rather saltier water from the Makassar Current, the extension of the Pacific North Equatorial and the Mindanao Currents, into the Java Sea, while fresher water from the South China Sea is transported eastward into the Java Sea during boreal winter monsoon season. This leads to one SSS peak usually in September and one minimum in January/February each year, which can be observed also in the time series for TA and DIC. The SST seasonality imprints its variability onto the time series of temperature dependent parameters such as solubility and the calcification process and, as a secondary effect, onto the others as well. Furthermore, in these Java Sea time series, the previously mentioned El Niño events are also visible (see peaks in SST, calcification rate and pCO 2 , more negative CO 2 flux (more outgassing) and low pH values in 1973, 1983, 1997/1998).
Simulated calcification rates are shown in Table 2 It is obvious from Figures 7, 10 that the average net production of aragonite per unit area is highest in the Gulf of Thailand followed by the one in the Java Sea. However, due the smaller area, the total regional amount of calcified aragonite in the Gulf of Thailand is lower than in the Java Sea. The lowest net production occurred in the Malacca Strait falling below the one of the entire Sunda Shelf.
The carbonate precipitation did not effect the long-term trends described above but increased the pCO 2 and thus the CO 2 outgassing as expected. Nevertheless, in 2007, the increasing atmospheric CO 2 concentration exceeded for the first time the oceanic and turned the Java Sea into a CO 2 sink for a couple of months.
FIGURE 7 | Simulated parameters as spatially averaged yearly running means for the entire Sunda Shelf Sea (gray), Gulf of Thailand (red), Malacca Strait (green), and the Java Sea (blue). pH is in seawater scale; "Omega" is the saturation state of aragonite; "calc rate" means aragonite calcification rate. Locations of regions are shown in Figure 2C.

Carbonate System at Certain Points
For spatially non-averaged simulation results, seven different locations have been selected to present the differences of the simulation results "without-with" calcification, again as yearly running means, see Figure 9: two neighboring points (nBan) in the Banten Bay, just north of west Java Island, with 55 and 10% coral reef coverage, two almost neighboring points (wMak) in the Spermonde Archipelago, southern Makassar Strait, with 44 and FIGURE 8 | Simulated parameters as spatially averaged monthly means for the Java Sea without biological processes (blue) and with calcification (red). pH is in seawater scale; "Omega" is the saturation state of aragonite; "calc rate" means aragonite calcification rate. SSS and SST are for both simulations.
14%, and three points without any coral reefs: central western Java Sea (wJaS), central and southern Malacca Strait (cMlc, sMlc). All locations are marked in Figure 1.
The time series for both nBan locations nearly match each other, even though their reef coverage with 55 and 10% is quite different. Their calcification rates per m 2 are nearly the same because of similar water properties at both locations-due to the current system and the resulting dispersion. The calcified mass (not shown), however, is much larger at nBan55 due to the higher coral occurrence. In contrast to this, both wMak points (44 and FIGURE 9 | Differences of simulated parameters without-with reef calcification as yearly running means for certain grid points (see Figure 1). Colors are mapped to locations in the calc rate graph: nBan: north of Banten Bay, wMak: west of City of Makassar, wJaS: western Java Sea, cMlc/sMlc: central/southern Malacca Strait, numbers in location names mean % of reef coverage. pH is in seawater scale; "calc rate" is aragonite calcification rate (only simulation with corals). CO 2 flux and pH values are shown twice: the lower two panels are monthly values to visualize the seasonal variability.
14%, green and blue curves in Figure 9), which are not direct neighbors, show slightly different but comparable calcification rates per m 2 , but their unequal calcified masses, due to different coral occurrence, influence the water properties in their grid cells also quite differently. As expected, the changes in the wMak44 water are larger for all parameters, e.g., with higher pH decrease of around 0.025 at wMak44 vs. 0.01 at wMak14, higher pCO 2 increase of 15-25 µatm at wMak44 vs. 5-15 µatm at wMak14. This leads also to a stronger increased CO 2 release at wMak44 than at wMak14. Naturally, dispersion of properties plays also an important role in marine environments. This is evident from the time series of the three locations without any coral coverage (wJaS00, cMlc00, sMlc00), where the visible differences between simulation results without and with calcification in Figure 9 were not locally generated but transported to these places.
To stronger illuminate the variability, the lower two panels of Figure 9 display monthly averages instead of yearly running means for two important parameters, the atmosphere-ocean CO 2 flux and the pH value. Seasonally and locally (e.g., at wMak44, green curves), the pH value can decrease by up to 0.06 units, and the CO 2 outgassing can increase by up to almost 10 mmol C m −2 d −1 .

DISCUSSION
In this section, we discuss our model results according to the most important processes we aimed to investigate in this first step of our model application: the air-sea flux of CO 2 , the ocean acidification, and the coral reef calcification. They are connected to each other, because the first one has implications on the second and subsequently on the third one.
In the previous section, we have shown that our model results agree quite well with observational data already for the first part of simulation, which was without any biological processes. This might suggest that the basic biological processes in the tropical shelf sea-photosynthesis, respiration, calcification, dissolution-are overall well balanced. The consideration of realistic net coral reef calcification rates leads to a slight change of geochemical parameters in the seawater and in the exchange of CO 2 between ocean and atmosphere. Obviously, this has an overall moderate impact on the carbonate system in the tropical shelf sea, which might, however, alter with climate change effecting the several involved processes differently.

Air-Sea Fluxes of CO 2
For the entire Sunda Shelf, the air-sea CO 2 exchange are generally directed from ocean to atmosphere with a seasonal variation due to the seasonality of the atmospheric CO 2 concentration. In some areas like the Gulf of Thailand, emission of CO 2 occurs only seasonally, in other parts like the Java Sea, seasonal CO 2 uptake started in recent years. Obviously, there is a trend toward less outgassing or more uptake for all areas superimposed by the trend of growing atmospheric CO 2 content. Assuming constant trends in atmospheric CO 2 concentration, in SSS and SST development, according to Table 2, there will be no CO 2 emission in the Sunda Shelf domain any more after 30-35 yearswith all consequences on the ocean acidification and the marine biology.
The effect of Equation 3 showing the production of CaCO 3 and additional CO 2 in the seawater is a shift of the curve of airsea CO 2 flux toward more emission. For the Sunda Shelf water, this effect amounts to about 10% from 15 to 17 Tg C y −1 because of the overall small coral reef coverage. For local areas, however, the effect can be as high as visible in Figure 9 by the green curve showing an increase of outgassing by up to 3.5 mmol C m −2 d −1 for a location with a 44% coral reef coverage.

Ocean Acidification
The simulated trend for the pH values (see Table 2: decrease of around 0.1 per century) is consistent with publications for other regions (e.g., Santana-Casiano et al., 2007, Canary Islands, around 0.1 units per century). This indicates that biological processes may not change the trend significantly. Most likely, they influence the pH level, correspondingly to what we found in our results without and with coral reef calcification. Following DeCarlo et al. (2017), the consideration of biological processes in coral reefs would increase the pH value and also the aragonite saturation state because of net CO 2 fixation from net ecosystem production (NEP). This would, in turn, enhance NEC. But NEC would reduce the pH value again. Also Silverman et al. (2007) found from their observations a dependency of NEC on nutrient concentrations, which corresponds to a dependency of NEC on NEP.
Therefore, the consideration of biological processes other than reef calcification by means of NEP would, on the one hand, lead to an increase of the pH level and a subsequent increase of calcification rates, which would then, on the other hand, lower Regions are shown in Figure 2C; "inner Sunda Shelf" corresponds to the sum of all four colored regions in that figure.  (2010)(2011)(2012)(2013)(2014) and trends (1965-2014, 50 years, in units per century) in simulated time series without biology ("no bio") and with calcification ("w calc"). "i. Sunda Shelf" means inner Sunda Shelf consisting of all color-shaded sub-domains on the Sunda Shelf as presented in Figure 2C. "flCO 2 " is the air-sea CO 2 flux (neg.: outgassing), "ratea" the aragonite calcification rate, "cf" calcifying fluid. *These aragonite calcification rates are the sum for the entire Sunda Shelf as depicted by the 200 m line in Figure 2C, where the other region locations are also shown.
the pH level again. The estimations of the net effect of NEP plus NEC, i.e., the balanced state at the end including the rates of NEP and NEC, was not possible within the current study. However, it would definitely be a further research project worth to be performed. In other intensively studied shelf areas like the North Sea, biological processes such as CO 2 uptake during photosynthesis and its subsequent export into sub-surface water or directly into the underlying sediment strongly influence the air-sea fluxes of CO 2 (Lorkowski et al., 2012). In their studies, this function was declining during 1970-2006. Others like Metzl et al. (2010) found also declining carbon shelf pumping effected by increasing SST and decreasing pH values. Similar results were obtained by Lauderdale et al. (2016). These authors quantified the air-sea CO 2 flux drivers and find an oceanic CO 2 uptake due to biological activity, opposing all other drivers, for the tropical Indian, Pacific and Atlantic Oceans. However, their findings result from investigations in the open ocean and are not necessarily valid for a tropical shelf sea like the Sunda Shelf region.

Coral Calcification
Our estimated calcification rate for the entire Sunda Shelf amounts to 5.5 Tg C y −1 ( Table 2). Bauer et al. (2013) give a rate of 150 Tg C y −1 for the global coastal inorganic carbon burial, Chen and Borges (2009) publish a rate of 180 Tg C y −1 for the global particulate inorganic carbon burial including the fish catches on shelves. Therefore, according to our model results, the contribution of the coral reef calcification only on the Sunda Shelf is 3-4% of the global coastal carbon extraction. If we take this percentage for all warm water coral reefs and assume that the reef area on the entire Sunda Shelf amounts to 10% of the global coral area (see Introduction and Table 1: 43,128 out of 423,582 km 2 ), the global contribution of warm water coral reef calcification to the permanent extraction of inorganic carbon from the water column would be 30-40%.
For the simulation of the calcification process in a coral reef, we applied an equation for the net calcification rate, which represents the gross calcification minus biological and chemical dissolution (NEC of the entire coral reef communities). However, dissolution of CaCO 3 is not necessarily a certain and constant proportion of the calcified aragonite mass as published for example by Silverman et al. (2007). Therefore, an improvement of the model results would be a consideration of dissolution rates to make this process depend on the aragonite saturation and the temperature in the seawater, possibly also on the nutrient concentration.
As obtained from the model results, the trend of SST is an increase by more than 1 • C per 100 years, which has an enhancing effect on the calcification process (Cooper et al., 2012). The pH values are decreasing by 0.07-0.11 units per 100 years, which is lowering the calcification rate (Kleypas et al., 2006;Hoegh-Guldberg et al., 2007). Hence, these two trends have an opposing effect on the calcification process. However, according to our model results, the dominant component is the pH value reducing the availability of CO 2− 3 in the calcifying fluid when decreasing in the seawater. Additionally, the SSS time series show also a decreasing trend (1-2.8 psu per 100 years), due to increasing precipitation (rain), reducing the availability of Calcium necessary for calcification.

SUMMARY AND CONCLUSIONS
A numerical model system was applied to the tropical Indonesian Seas to simulate the CO 2 air-sea flux under the influence of increasing atmospheric CO 2 and under the additional influence of aragonite calcification from coral reefs for a period of almost 60 yeas. As a first step, simulations were performed without consideration of any biological processes. The long-term time series show realistic values and trends for relevant geochemical parameters of the physical and carbonate systems, supporting the findings of other authors regarding ocean acidification. The good agreement of these simulation results (no biological processes) with estimations from observational data suggests that biological processes taking up and releasing CO 2 are currently well balanced in these tropical regions. As a second step, we repeated the simulation but considered reef calcification, achieving realistic results for the calcification rates (around 1,400 g CaCO 3 m −2 y −1 ). The effect of this was an increase of CO 2 emissions. During the period from 2010 to 2014, calcification enhanced the average CO 2 emission of the Sunda Shelf by more than 10% from 15 to 17 Tg C yr −1 due to lowering the pH and increasing the pCO 2 in surface water. The long-term simulation results for the total rate of reef calcification as well as a lowering of reef calcification rate.
The Sunda Shelf region used to be generally a source of CO 2 for the atmosphere. However, influenced by increasing sea surface temperatures, enhancing outgassing, and increasing atmospheric CO 2 concentrations, damping outgassing, the CO 2 emissions were reduced on the Sunda Shelf. In some parts, the CO 2 flux was even reversed to oceanic uptake. Overall, the entire Sunda Shelf is expected to turn into a permanent atmospheric CO 2 sink within the next 30-35 years, assuming that our estimated trends are constant. On the long-term, this definitely would have a severe impact on the marine biological environment.
One of our major findings from our simulations is that tropical marine environments are really not all the same. It is possible and necessary to distinguish different regimes even within the limited Sunda Shelf region. They show different influences from calcification, different behavior in the carbonate system with some regions acting as a CO 2 source and others as a CO 2 sink for the atmosphere, and-more important-different trends in the long-term time series. This points out that it is important but not sufficient to investigate the global impact of climate change, because consequences on the regional scale might significantly defer from the global overall results and be much more severe.
To obtain further realistic and reliable results, especially for climate projections, the simulation of net reef calcification should consider both aragonite production and dissolution rates, not only an overall net rate, as well as primary production and respiration processes, because the several components might change differently with climate change. Furthermore, a temperature dependent limitation of coral growth will be incorporated to enable a consideration of further relevant processes like coral bleaching. This, however, depends on the availability of a comprehensive data inventory for this region, which provides corresponding and consistent biogeochemical and nutrient data to initialize and force the numerical model.

AUTHOR CONTRIBUTIONS
BM: Did most of research published. TR: Gave worthy research questions and expertise recommendations, took part in writing the manuscript. TP: Gave worthy research questions and expertise recommendations, did manuscript reviews, PI for the project.

ACKNOWLEDGMENTS
This project was funded by the German Ministry of Education and Research under the grant 03F0642D in the frame of the German-Indonesian cooperation SPICE III.
We would like to thank the Japan Meteorological Agency for hosting the website of World Meteorological Organization's (WMO) World Data Centre for Greenhouse Gases (WDCGG) and also the National Oceanic and Atmospheric Administration (NOAA) of the U.S.A. for providing the greenhouse gas data (measured atmospheric CO 2 ) via http://ds.data.jma.go. jp/gmd/wdcgg or ftp://aftp.cmdl.noaa.gov/data/greenhouse_ gases/co2/flask. Their data were basis to derive our formula for the atmospheric CO 2 concentration in our model simulations.
Furthermore, we would like to thank the Carbon Dioxide Information Analysis Center (CDIAC), located at the U.S. Department of Energy's (DOE) Oak Ridge National Laboratory (ORNL), for providing global oceanic carbonate data (TA, DIC) via their website http://cdiac.ornl.gov.
We are also thankful to Francisca Wit (Leibniz Centre for Tropical Marine Research, Bremen, Germany) for offering their measured data collected during our joint research cruises in Oct 2012 and April 2013, which we used for comparison with model results.