Local oOceanic CO 2 outgassing triggered by terrestrial organic carbon ﬂuxes during deglacial ﬂooding

. Oceans play a major role on the e Exchanges of carbon with between the ocean and the atmosphere are key processes that inﬂuence and thereby on past climates with glacial/interglacial via glacial/interglacial variations of the CO 2 concentration. The melting of ice sheets during deglaciations lets the induces a sea level rise which leads to the ﬂooding of coastal land 5 areas, resulting in the transfer of terrestrial organic matter to the ocean. However, the consequences of such ﬂuxes on the ocean biogeochemical cycle and uptake/release of CO 2 are poorly constrained. Moreover, this potentially important exchange of carbon at the land-sea interface is not represented in most Earth System Models. We present here the implementation of terrestrial organic matter ﬂuxes into the ocean at the transiently changing land-sea interface in the Max Planck Institute for Meteorology Earth System Model (MPI-ESM) and investigate their effect on the biogeochemistry during the last deglaciation. 10 Our results show that during the deglaciation, most of the terrestrial organic matter inputs to the ocean occurs during Meltwater Pulse 1a (between 15-14 ka) which leads to the transfer of additional 21.2 GtC of terrestrial carbon origin (mostly originating from wood and humus) to the ocean. Although this additional organic matter input is relatively small in comparison to the global ocean inventory (0.06 %) and thus doesn’t have an impact on the global CO 2 ﬂux, the terrestrial organic matter ﬂuxes initiate oceanic outgassing at regional hotspots like in Indonesia for a few hundred years. Finally, sensitivity experiments highlight 15 that terrestrial organic matter ﬂuxes are the drivers of oceanic outgassing in ﬂooded coastal regions during Meltwater Pulse 1a. Furthermore, the magnitude of outgassing is rather insensitive to higher carbon to nutrients ratios of the terrestrial organic matter. Our results provide a ﬁrst estimate of the importance of terrestrial organic matter ﬂuxes in a transient deglaciation simulation. Moreover, our model development is an important step towards a fully coupled carbon cycle in an Earth System Model applicable for simulations atof glacial/interglacial cycles.

. Recent modelling estimates also indicate that the land carbon content was lower of 450 to 1250 GtC during the LGM than during the Preindustrial (e.g. Menviel et al., 2011;Jeltsch-Thömmes et al., 2019;Müller and Joos, 2020) due to colder and drier climate conditions and to the presence of large ice sheets. However, the exact mechanisms responsible for the glacial/interglacial atmospheric CO 2 changes remain unclear and all the previous modelling attempts to explain these variations only considered the ocean-atmosphere CO 2 gas exchange and carbon fluxes between the land carbon pool and the atmosphere, Stanford et al., 2006). One of this event, the Meltwater Pulse 1a (MWP1a) around 14.65 ka, is associated with a rapid sea level increase from 8.6 to 20.2 m within 500 years (Deschamps et al., 2012;Liu et al., 2016;Lin et al., 2021). The source of the large amounts of meltwater is, yet, unclear. Either a partial melting of the Northern Hemisphere Ice Sheets or of the Antarctic Ice Sheet, or changes of both ice sheets are currently investigated by modelling studies (Mackintosh et al., 2011;Golledge et al., 2014;Weber et al., 2014;Gomez et al., 2015;Gregoire et al., 2016;Yeung et al., 2019). However, due to the rapid sea level 75 change, MWP1a is a particularly interesting period to look at the role of terrestrial organic matter fluxes on atmospheric CO 2 variations. In this study, the freshwater inputs of MWP1a to the ocean are deduced from a prescribed ice sheet reconstruction and are located in the Northern Hemisphere (Tarasov et al., 2012). We are aware of the uncertainty of the origin of this meltwater pulse which could potentially impact the location of flooded land coastal areas, and ultimately the land-sea fluxes. 80 There is a lack of global coupled climate models that consider fluxes of terrestrial organic matter induced by flooding, as well as their consequences on the ocean over long time period such as deglaciations. Only few previous studies conducted time-slice sensitivity experiments with EMICs to quantify the impact of changes in marine physical and biogeochemical conditions on ocean outgassing. Running transient simulations with General Circulation Models remain challenging. We present here for the first time a coupled transient simulation over the last deglaciation using MPI-ESM (Max Planck Institute Earth System Model) 85 that takes into account (1) a fully interactive adaptation of the ocean bathymetry with corresponding changes of the land-sea distribution (Meccia and Mikolajewicz, 2018), (2) a transient river routing (Riddick et al., 2018) and the representation of new processes that are (3) the automatic adjustment of marine biogeochemical tracers under changing ocean bathymetry and land-sea distribution and (4) the fluxes of terrestrial organic matter between land and ocean during flooding events. A special focus is placed on the time period of the large freshwater inputs between 15-14 ka and corresponding to Meltwater Pulse 1a 90 (MWP1a) to address the role of terrestrial organic matter fluxes on the CO 2 fluxes between the ocean and the atmosphere in context of this millennial event.
2 Model description MPI-ESM is a comprehensive Earth system model composed of four components and one coupler: the atmospheric component ECHAM6.3, the ocean dynamics component MPIOM1.6, the ocean biogeochemistry component HAMOCC6, the land, 95 hydrology and dynamic vegetation component JSBACH3.2 and the coupler OASIS3-MCT. It is currently used in the version 1.2 following the latest updates detailed in Mauritsen et al. (2019). Several components of MPI-ESM have been further extended with new developments to take into account a varying land-sea mask during the period of sea level changes (Meccia and Mikolajewicz, 2018) and the river routing (Riddick et al., 2018), as described in the following sections. In addition, we present here new model developments that allow it to take into account the terrestrial organic matter fluxes between land and 100 ocean at a transiently changing land-sea interface as well as their effect on the ocean biogeochemistry.

HAMOCC
The global ocean biogeochemical model HAMOCC (HAMburg Ocean Carbon Cycle) is part of the ocean component MPIOM (Max Planck Institute Ocean Model) which is used to simulate the oceanic physics by resolving primitive equations under the hydrostatic and Boussinesq approximation on a C-grid with a free surface for every time step (1h) (Jungclaus et al., 2013). 105 HAMOCC simulates the biocheochemical tracers and processes in the oceanic water column, the sediment and at the air-sea interface. HAMOCC has been previously described in Ilyina et al. (2013) and then revised in Mauritsen et al. (2019) with additional processes like marine nitrogen fixation by Cyanobacteria (Paulsen et al., 2017) and nitrogen deposition. Recently, additional ocean tracers have been added in HAMOCC to simulate carbon isotopes (Liu et al., 2021). In this paper we use the HAMOCC6 version as described in Mauritsen et al. (2019) with specific new features for the coupling of the land-sea-110 continuum that are described in the next sections. The grid configuration of HAMOCC is identical to the ocean component MPIOM and consists of a bipolar grid with one pole over Greenland and another pole over Antarctica. The coarse-resolution ocean grid used for paleoclimate purpose is noted GR30 and has a horizontal resolution of about 300 km. The vertical resolution is composed of 40 unevenly spaced layers with level thickness increasing with depth. Details on the coupling between MPIOM and HAMOCC can be found in Maier-Reimer et al. (2005). dependent remineralisation rate with a Q10 of 2.3. The remineralization of organic matter can be aerobic or anaerobic with denitrification or sulfate reduction when O 2 concentration passes a critical threshold. This version of HAMOCC also uses an updated procedure to account for weathering fluxes entering the ocean. Instead of a globally uniform distribution, weathering fluxes are only distributed along the coast. The water column interacts with atmospheric components by air-sea exchange of gaseous tracers (O 2 , N 2 , N 2 O and CO 2 ). All atmospheric concentrations are prescribed. Dust deposition from the atmosphere 130 to the ocean is also taken into account based on the time slice estimates of Albani et al. (2016) (see details in Section 2.4).
The simulation of the oceanic sediment follows the approach of Heinze and Maier-Reimer (1999). The biogeochemical tracers of the sediment include dissolved inorganic carbon(DIC), alkalinity as well as dissolved phosphate, nitrate, silicate, oxygen, iron, N 2 and H 2 S for the pore water fraction and detritus, opal, calcium carbonate and clay for the solid fraction. The processes 135 of decomposition of detritus, dissolution of opal and calcium carbonate and the carbon chemistry are similar to those in the water column. There is vertical diffusion of dissolved tracers in pore water within the sediment column, as well as a diffusive exchange with the water column above. The sediment is resolved by 12 layers with increasing thickness and decreasing porosity from top to bottom, the bottom layer being defined as an underlying burial layer. In the burial, no biogeochemical processes are active and there is no pore water diffusion between burial and active sediment layers. This HAMOCC version 140 also includes erosion under higher bottom shear stress. At most half of the first sediment layer volume of the solid components can be eroded within one time step. Due to remineralization and dissolution of solid components as well as erosion, gaps might be produced within the sediment layers. They are removed by shifting the solid sediment constituents downward within the sediment column. To guarantee that each sediment layer is properly filled with solid components according to solid fraction, the sediment column is shifted upward and filled with sand from the burial layer. Sand is an inert component and can not be 145 eroded (Mathis et al., 2019).

Dynamic land-sea mask and hydrological discharge
In context of ice sheets growth or decay and associated changes in freshwater input and sea level, it is necessary to take into account the variations in ocean bathymetry and coastlines when performing transient simulations. Meccia and Mikolajewicz (2018) describe a new automatic method for interactive bathymetry and land-sea mask changes in MPIOM, allowing simula-150 tions for transient time periods. To do so, an ice sheet reconstruction is used to compute the time-dependent freshwater fluxes to the ocean and to derive the topography, which is then used to obtain changes in bathymetry and land-sea mask. Changes in ocean bathymetry and the land-sea mask are updated every 10 years. In the case of ocean expansion, new grid cells are flooded by the surrounding water. For HAMOCC, this also includes pore water of new sediment grid cells. The solid fraction of the new sediment column is filled with sand. This automatic bathymetry adjustment conserves mass (i.e. salt content and all other 155 oceanic tracers) at the global and regional scales (changes in ocean volume match the freshwater input and the global inventory of tracers is constant in the absence of sources or sinks).
Changes in ice sheets volume since the Last Glacial Maximum (LGM) induce changes in the river pathways that need to be taken into account to properly estimate the runoff during the last deglaciation. Riddick et al. (2018) presented (a) new developments within the land surface component JSBACH (Jena Scheme for Biosphere-Atmosphere Coupling in Hamburg) (Reick et al., 2013;Kleinen et al., 2020;Reick et al., 2021) and the Hydrological Discharge model (Hagemann and Dümenil, 1997;Hagemann and Gates, 2001) to automatically generate dynamic river directions and flow parameters for past conditions. First the orography is adjusted and corrected. Then river directions and flow parameters are generated and the location of the river mouths are determined. The ice sheet height and isostatic adjustments are taken from an ice sheet reconstruction and 165 the land-sea mask is generated using the technique described above. These changes in freshwater discharge to the ocean have significant impacts on the global oceanic circulation (via the North Atlantic and Arctic Oceans) and as a consequence also affect atmospheric circulation.

Land-sea carbon and nutrients transfer
In this version of MPI-ESM, we account for the fact that coastal areas are flooded due to increasing sea level (following the 170 melting of the ice sheets) (coastal areas are flooded). ECHAM6.3 does not include fractional grid boxes, i.e. a single grid box is treated as land (ocean) as long as the land (ocean) fraction is larger than 50 %. Flooding, i.e. conversion from land to ocean, occurs if the land fraction within one grid cell is less than 50 %. Drying of an ocean grid cell in the case of a decrease of the sea level is also considered. In this case, all pore water tracers, i.e. phosphate, nitrate, dissolved inorganic carbon and alkalinity are redistributed to the water column to guarantee mass conservation and the solid parts of the sediment are considered as inactive 175 land pools.
In the land component JSBACH, each grid box represents the diversity of plant functional types (PFTs) depicting various vegetation types (Reick et al., 2021). 11 PFTs are described in the model: 8 for natural vegetation and 3 for land-use types (which are not used in this deglaciation simulation). The land module also accounts for coupled carbon and hydrological cycles. As part of the atmospheric component, the land component JSBACH uses the same coarse-resolution grid noted T31 with an 180 approximate grid spacing of 400 km (different from the GR30 grid of MPIOM). In the case of flooding, the terrestrial organic matter of the flooded grid cell is collected on the T31 grid and remapped to the GR30 grid. The transfer from land to ocean is mass conserving.
As a consequence of flooding events, the land carbon pools (wood, woody litter above and below ground, humus) are trans-185 ferred to corresponding oceanic pools in the water column and sediment ( Figure 1). So in addition to the 25 prognostic tracers already existing in HAMOCC, we added 3 new prognostic tracers to represent wood, woody litter and humus. In case of flooding, woody litter above ground is distributed uniformly over the water column of the corresponding ocean grid cells. Woody litter below ground and humus are transferred to the sediment where they participate in all sedimentary processes and can be eroded ( Figure 1). Wood is located at the water-sediment interface. It is not advected and does not participate in erosion. All 190 terrestrial organic particles in the water column sinks with the same speed as marine particulate organic matter. All terrestrial organic matter is remineralized using oxygen in case of aerobic remineralization or NO 3 and NO 2 if anaerobic remineralization takes place (see equations below). The prescribed remineralization timescales in water are 100 years for wood, 10 years for woody litter and 5 years for humus. It is crucial to note that the stoichiometry of terrestrial organic matter differs from that of marine organic matter defined as C:N:P = 122:16:1 (Takahashi et al., 1985). The carbon to nutrients ratio for woody litter 195 (above and below ground) is defined as C:N:P = 7600:51:1 (Goll et al., 2012). The carbon to nutrients ratio for wood is set to C:N:P = 3650:11:1 (Goll et al., 2012). The carbon to nutrients ratio for humus is C:N:P = 465:10:1 (Goll et al., 2012). Remineralization of terrestrial organic matter leads to changes in dissolved inorganic carbon, phosphate, nitrate, dissolved oxygen and alkalinity in the water column or the sediment. In contrast to the long-living terrestrial material that is transferred to the ocean, short-living material from the vegetation (green biomass, non-woody litter above and below ground) is recycled within 200 one year. We assume conversion of land to ocean by flooding happens over 10 years and after this time, all the short-living terrestrial organic matter is remineralized and emitted as CO 2 to the atmosphere (Figure 1). considered to be remineralized during the flooding event and is treated as a CO 2 source for the atmosphere (Figure 1). However, since we used prescribed CO 2 concentrations in this simulation, the flooding induced terrestrial carbon emitted to the atmosphere from short-living land sources has no effect on the climate. Assuming that terrestrial carbon compounds that are carbohydrates (CH 2 O) x give the classical organic matter composition defined as (CH 2 O) x (NH + 3 ) y H 3 PO 4 , we use the equations of Paulmier et al. (2009) presented below to calculate the oxygen and nitrate demand for remineralization, as well as alkalinity change. Considering that land biomass contains no excess H + compared to organic material in the ocean we obtain: for oxygen We use the equation for aerobic remineralization with an autotrophic shortcut to nitrate as defined in Paulmier et al. (2009): The change in alkalinity is given by (-d-1) following Eq. (2). For complete denitrification (with conversion of the ammonium 215 produced by anaerobic remineralization into N 2 ) we use the following equation (Paulmier et al., 2009): For the NO − 3 change, the value corresponds to (alkalinity change + 1), i.e. ( 4 5 a + 1 5 b -2 5 c + 1). This leads to the overall values presented in Table 1 for the different carbon and nutrients compositions of terrestrial organic matter.

Model and experiment setup 220
Two spin-up runs (between 26-24 ka and 24-21 ka) were first performed to bring the physical, climatic and biogeochemical systems into quasi-equilibrium. We then ran the deglaciation simulation starting at 21 ka. The initial conditions and forcing are described next. The greenhouse gases concentrations, i.e. CO 2 , CH 4 and NO 2 are from Köhler et al. (2017) and the orbital parameters are taken from Berger and Loutre (1991). The atmospheric CO 2 concentration is prescribed during these simulations to test and validate the state of the model with the new developments presented in previous section. A new simulation with 225 prognostic atmospheric CO 2 will be performed in the future after additional model development to address the gap on the interaction between the ocean biogeochemistry and the climate during the last deglaciation. The model is also forced with the ice sheet reconstruction GLAC-1D (Tarasov et al., 2012) and dust deposition from Albani et al. (2016) which is linearly interpolated over the deglaciation since the data is only available for specific time-slices (21 ka, 16 ka, 10 ka, 8 ka, 6 ka, 4 ka, 2 ka). The initial marine nutrients and carbon inventories are set to similar values to those of the present -day ocean with global surface mean concentrations of 4.2 mmol N m −3 for nitrate, 0.52 mmol P m −3 for phosphate, 9.2 mmol Si m −3 for silicate and 1915 mmol C m −3 for DIC. We didn't adjust the nutrient concentrations for the 3.5 % change in the oceanic volume between the present day and the LGM. The global mean alkalinity concentration is set higher than present-day with 2300 mmol m −3 to stabilize the marine chemistry to the lower LGM atmospheric CO 2 concentration. Weathering fluxes are set to a global constant value for the entire simulation to compensate for burial fluxes, based on the previous spin-up runs. The global weathering rates 235 used to calculate the coastal input fluxes correspond to 1 kmol P s −1 for organic material, 1720 kmol C s −1 for inorganic material and 210 kmol Si s −1 for silicate. The weathering fluxes are independent from the river discharge even if they might locally show small variations because the length of the coastline varies. Moreover, we consider that high latitudes coastlines get only 20 % of the global value to avoid excessive inputs of silicate that would occur if the weathering rates were distributed uniformly over the coastlines (Lacroix et al., 2020). Finally, we assume that sea-ice contains only pure water and thus does not 240 carry any biogeochemical tracer or salt. The land-sea mask, bathymetry, river directions and flow parameters are updated every 10 years.
Different sensitivity experiments have also been performed during Meltwater Pulse 1a by branching off new simulations at 15 ka from the deglaciation simulation. These sensitivity experiments ran for 1000 years. One sensitivity experiment has been run without terrestrial organic matter fluxes to the ocean to investigate the changes in sea-air carbon flux. Furthermore, a A 245 set of two sensitivity experiments has been run with modified carbon to nutrients ratios for terrestrial organic matter entering the ocean (details in Section 3.3). Furthermore, one sensitivity experiment has been run with modified stoichiometry and remineralization rates for terrestrial organic matter. reconstruction. In our model, the ice sheet volume decrease is considered as liquid meltwater discharge, ignoring the discharge 285 of icebergs which would lead to slower freshwater input to the ocean. Thus, a pulse-like meltwater occurs during 15-14 ka, leading to a rapid AMOC weakening. The variability of the simulated AMOC is mainly regulated by temporal changes in the volume of the prescribed ice sheet reconstruction. In our model, the ice sheet volume decrease is considered as liquid meltwater discharge, ignoring the discharge of icebergs. Freshwater inputs deduced from the GLAC-1D ice sheet reconstruction show low variations during the LGM and only a slight increase during the Heinrich Stadial 1. Thus, we can't expect pronounced AMOC 290 changes during the period of Heinrich Stadial 1 in the model. Between 15-14 ka, we simulate a decrease of the AMOC strength following massive freshwater inputs in the Northern Hemisphere. This period of MWP1a is also characterized by a rapid sea level increase, which is recorded in radiocarbon dates from the Sunda Shelf and U/Th measurements on corals offshore from Tahiti, confirming a timing of MWP1a between 14.65 to 14.31 ka (Hanebuth et al., 2000;Deschamps et al., 2012). However, the temporal variation of the AMOC strength estimated from 231 Pa/ 230 Th tends to show already a decrease between 17.5-15 295 ka, i.e. before the MWP1a, and an increase back to a high value between 15-14 ka (McManus et al., 2004). To achieve a good agreement between simulated and proxy-data derived AMOC variations, He (2011) showed in a modelling study the necessity of meltwater fluxes from Antarctica during MWP1a. Other hosing experiments also emphasize the sensitivity of the oceanic circulation, and thus the AMOC, on the location of the freshwater input (e.g. Smith and Gregory, 2009;Menviel et al., 2011).
As previously discussed, the meltwater input deduced from GLAD-1D is located primarily in the Northern Hemisphere, which 300 might explain the different temporal evolution of the simulated AMOC. In the following, we refer to MWP1a as the period of rapid sea level change due to large freshwater inputs to the ocean to evaluate the effect of land-sea exchanges during flooding events.
The sea level change over the last deglaciation can be estimated based on the freshwater inputs to the global ocean, that are set accordingly to the GLAC-1D ice sheet reconstruction. Considering that the sea level at 21 ka is similar to the estimated 305 one from Spratt and Lisiecki (2016)   For the land part, the total carbon in all terrestrial carbon pools (litter, soil and vegetation) increases from 922.9 GtC to 1302.7 315 GtC between 21-15 ka and reaches 1563.6 GtC at 12 ka (Figure 2eg). This increase is explained by the progressively warmer climate during the deglaciation and by the increase of C 3 plants and Gross Primary Productivity (GPP) on land. Even if our simulation currently doesn't continue beyond 12 ka, the total land carbon content evolution since the Last Glacial Maximum is within the range of other modelling studies like Prentice et al. (2011), O'ishi and Abe-Ouchi (2013), Ganopolski and Brovkin (2017), Jeltsch-Thömmes et al. (2019) and Müller and Joos (2020). Furthermore, in a setup with a fully coupled land and ocean 320 carbon cycle it would be the ocean that supplies this additional land carbon, which would demand a constant outgassing over the period during which the land carbon inventory increases. During MWP1a, the C 3 cover fraction and the land GPP decrease, indicating the replacement of trees by grasslands and resulting in a decrease of 62.38 GtC stored in litter and vegetation pools as shown in Fig. 2eg. 325 We can also evaluate the ability of the model to reproduce the biome distribution since the LGM before the beginning of MWP1a at 15 ka in comparison to pollen data like the BIOME6000 Version 4.2 reconstruction (Harrison, 2017) based on the Palaeovegetation Mapping Project (Prentice and Webb III, 1998;Prentice et al., 2000;Harrison et al., 2001;Bigelow et al., 2003;Pickett et al., 2004). To do so we used the biomisation technique developed presented in Dallmeyer et al. (2019) to convert the different PFT cover fractions modelled in JSBACH into 9 biomes. We also used the best neighbour score (BNS) 330 metric method presented in Dallmeyer et al. (2019) to quantify the similarity between the modelled biomes and the pollen data from the BIOME6000 database. This method uses the surrounding grid boxes of the studied grid cell (in each direction of the T31 grid) to compare with the pollen record. The agreement for each record is calculated with the distance weight of the best neighbour in each neighbourhood (using a Gaussian function) and varies between 1 if the modelled biomes in the grid box indicates the same biome as reconstructed and 0 if all grid cells in the neighbourhood disagree with the record. The BNS is 335 the mean of all individual neighbourhood scores. The LGM modelled biomes on Fig. 4a show an overall good agreement with the pollen data with a total BNS value of 0.52. At high latitudes of the Northern Hemisphere, tundra and boreal forests are simulated in regions that are not covered by ice, which is consistent with the few pollen datasets available at these locations (BNS value of 0.78 and 0.19 respectively). Temperate forest is modelled over part of North America, grassland over Europe and temperate/warm forest over East Asia. This is in agreement with the pollen record even if some local discrepancies are Although the LGM conditions were different from those at 15 ka before MWP1a, in absence of other global reconstructions we also used the LGM BIOME6000 pollen record to compare to model results. According to our model, the biome distribution 345 doesn't change much between 21 and 15 ka (Figure 4a,b) so that for many regions, the LGM pollen data show the same pattern as the simulated biomes at 15 ka. The BNS value at 15 ka is similar to the one at 21 ka for tropical forest, warm forest, savanna and desert. However, climatic differences between these two periods lead to small differences between the simulated biomes at 15 ka and the 21 ka pollen data which explains the lower total BNS value of 0.45 compared to 0.52 (Figure 4a,b). Part of the LGM tundra at high latitudes of the Northern Hemisphere is replaced by the boreal forest or grassland at 15 ka. At low latitudes, there is a slightly larger extent of the temperate forest over East Asia and of the tropical forest over South America at 15 ka. The tropical forest over Indonesia is however already present since the LGM. the BIOME6000 pollen record shows a good overall agreement with the modelled biome reconstruction (Figure 4). At high latitudes of the Northern Hemisphere, both tundra and boreal forests are modelled for the areas that are not covered by ice, which is consistent with the few pollen data sets available at these locations. Temperate forest is modelled over part of North America, grassland over Europe and    (Harrison, 2017). The right plots indicate the best neighbour score, i.e. the similarity between the modelled biomes and the pollen data, for both time period.

Land-sea carbon fluxes
The transient adaptation of the land-sea mask and bathymetry during the deglaciation allows for flooded coastal land areas to be accounted for when the sea level changes. During the deglaciation, before and after Meltwater Pulse 1a (i.e. 21-15 ka and 14-12 ka), flooding is observed on almost all continental coasts but with a higher number of flooded land coastal regions located above 60°N ( Figure 5). During MWP1a, local coastal land surfaces are flooded in East Asia, Indonesia and Australia, 375 but again with a higher number of the flooded grid cells in nNorthern Europe since the major source of the meltwater is in the Northern Hemisphere ( Figure 5). In terms of area, the flooded high latitudes coastal regions represent 1.03 x 10 6 km 2 , which is larger than the 8.04 x 10 5 km 2 of flooded coastal area in Indonesia.  Figure 6). There are also differences in the size of the carbon fluxes depending on the origin of the terrestrial organic matter.
During MWP1a, 7.7 GtC are emitted to the atmosphere (related to the green vegetation), which corresponds to 26.6 % of the total terrestrial organic matter inputs (Table 2). We however refrain from any discussion since these fluxes don't have an impact on the climate in this simulation with prescribed CO 2 concentrations. 21.2 GtC goes to the ocean, i.e. 73.4 % of the total terrestrial organic matter inputs (Table 2).  Looking at the origin of the terrestrial organic matter inputs to the ocean, wood and humus are the largest contributors with 31.1 and 29.1 % of the total inputs respectively (atmosphere and ocean) and have even higher contribution of 42.5 and 39.6 % for the ocean part only. In comparison to the entire deglaciation between 21-12 ka, the flooding induced terrestrial carbon emissions terrestrial organic matter inputs to the atmosphere and the terrestrial carbon inputs to the ocean represent a similar proportion than during MWP1a with respectively 28.4 and 71.6 % for a total of 88.3 GtC (Table 2). Wood and humus are again 390 the major contributors of these land-sea fluxes. We have to emphasize that these numbers are for one millennium and thus the total terrestrial organic matter contribution to the ocean of 21.2 GtC is non negligible compared to the entire deglaciation Among all contributors of terrestrial organic matter during MWP1a, the Indonesian region is the largest one with inputs to the water column that represent 66.4 % of the total inputs ( Figure 7a) and explained by wood inputs. The second largest contributor during this meltwater pulse is Australia with 17.3 %, followed by East Asia, the high latitudes of the Northern Hemisphere, South America and Europe with respectively 11.1, 2.2, 1.8 and 1.2 %. All these terrestrial inputs are dominated by wood inputs.

400
As for the part that goes to the sediment, Indonesia and East Asia regions are also major contributors with 43.7 and 18.7 % of the total terrestrial inputs to the sediment during MWP1a (Figure 7b), mostly due to humus. Other regions subject to flooding like Australia, the high latitudes of the Northern Hemisphere and Europe also contribute to terrestrial organic matter inputs to the sediment in smaller amounts. So even if most of the flooded areas are located in the Northern Hemisphere during MWP1a, their contributions of terrestrial inputs to the ocean are less significant than equatorial and low latitudes flooded areas. These results are however only representative for a short time period relative to the rapid change induced by the meltwater pulse and could be different from the rest of the deglaciation. To investigate this, we also present the contribution of these same areas during the full deglaciation (21-12 ka). In comparison to the specific case of MWP1a, Europe, South America, East Asia and the high latitudes of the Northern Hemisphere have similar contributions of terrestrial carbon and nutrients to the water column over the entire deglaciation (Figure 7a). In contrast, North America (East and West) and land areas in the 410 high latitudes of the Southern Hemisphere show organic matter inputs of a few percent during the deglaciation, but not during the large meltwater pulse. But the largest contribution of terrestrial organic matter to the water column observed during the deglaciation, i.e. Indonesia, occurs during the large meltwater pulse. In terms of regional contribution to the total terrestrial inputs to the sediment, Indonesia is once again the largest contributor with 43.7 % (Figure 7b). East Asia contributes with 18.7 %, similar to Australia, but three times more important than Europe and the high latitudes of the Northern Hemisphere 415 with 6.5 %. Among the smallest contributors are North and South America between 1 and 3 %. Once in the sediment, part of the terrestrial organic matter is buried, i.e. 0.042 and 0.168 GtC respectively for the woody litter and humus during the entire deglaciation. Inputs to the water column due to sediment erosion are of the order of 0.0026 and 0.026 GtC for the woody litter and humus respectively. The remaining part is remineralized.

Implications for the ocean biogeochemical cycle 420
To assess the role of the terrestrial organic matter inputs on ocean biogeochemistry during MWP1a, we performed a sensitivity experiment without terrestrial organic matter fluxes between land and ocean (as described in Section 2.4). The sensitivity experiment is branched off at 15 ka from the reference deglaciation run. Figure 8 shows anomalies of surface alkalinity, dissolved inorganic carbon and phosphate between simulations with and without terrestrial organic matter inputs (i.e. deglaciation reference run -sensitivity experiment) as average over MWP1a. An equatorial box between 15°N and 15°S around Indonesia is 425 defined and subdivided in three parts: northern part, central part and southern part, to discuss the largest differences (see Figure   9).     Figure 9. Same as Fig. 8a, b, c but with a focus over Indonesia for the surface alkalinity (a), surface dissolved inorganic carbon (b) and surface phosphate (c) between the two simulations with and without terrestrial organic matter fluxes averaged over MWP1a. An equatorial box is defined between 15°N and 15°S and subdivided in three areas: northern part, central part and southern part.
Due to the stoichiometry of the terrestrial organic matter being higher than of the marine organic matter (C:N:P = 3650:11:1 for the wood, 7600:51:1 for the woody litter and 465:10:1 for the humus in comparison to 122:16:1 for marine organic matter) there is an excess of carbon in the ocean once the terrestrial organic matter has been remineralized. Indeed, net primary production uptakes only 122 moles of carbon per mole of phosphate which leaves a large fraction of remineralized terrestrial carbon in the water column. This leads the ocean to outgas CO 2 to the atmosphere. Similarly to what has been presented above, we 460 investigate the effect of terrestrial organic matter inputs on the CO 2 fluxes between ocean and atmosphere during MWP1a. Figure 10 shows the anomaly between the two simulations (with and without land-sea fluxes) of the mean surface CO 2 flux between the ocean and atmosphere. Indonesia shows differences with higher CO 2 outgassing values of up to 2.2 x 10 −9 kg C m −2 s −1 for the simulation taking into account the land-sea fluxes. These larger values originate from the remineralization of the wood and humus associated with tropical forest developed at that time, as explained in Section 3.1 and shown in Fig. 4. organic matter fluxes, like in East Asia with 0.17-0.22 x 10 −9 kg C m −2 s −1 or in the nNorthern Australian coast with 0.36-0.39 x 10 −9 kg C m −2 s −1 . Overall, the terrestrial organic matter fluxes only have a local influence on the ocean behaviour To investigate the influence of the carbon to nutrients ratio of terrestrial organic matter on the CO 2 fluxes between the ocean and atmosphere during MWP1a we performed two additional sensitivity experiments following the procedure described in Section 2.4. The first sensitivity experiment uses a higher C:N:P ratio of terrestrial organic matter, arbitrarily defined as twice the carbon and nitrogen values of Goll et al. (2012) used in the reference simulation since no higher estimations exist in the literature. The second sensitivity experiment uses a lower C:N:P ratio of terrestrial organic matter based on studies from Lerman et al. (2004) and Cleveland and Liptzin (2007). These changes in the stoichiometry also require an adjustment of the O 2 and NO − 3 demand 475 and corresponding ∆Alk changes during remineralization of terrestrial organic matter. All values are summarized in Table 3. We compare the largest outgassing events over Indonesia for the different simulations with and without terrestrial organic matter as well as for the simulations with different carbon to nutrients ratios for this organic matter. As already observed in Fig.   10, the outgassing of the ocean in the equatorial region of Indonesiato the atmosphere from North to South Indonesia during MWP1a is driven by the terrestrial organic matter fluxes during flooding events. We observe a CO 2 flux to the atmosphere 480 happening concurrently to the flooding events and then slowly decreasing with time until all the terrestrial organic matter has been remineralized. In nNorthern partIndonesia (12.5°N -96.5°E), the outgassing starts at 14.64 ka, reaches 3.8 x 10 −9 kg C m −2 s −1 and then decreases for 200 years which corresponds to the decay time of the wood material (i.e. stems) in the ocean ( Figure 11a). This outgassing event is due to the transfer during the flooding event of terrestrial material to the ocean, dominated by the wood with 1.72 GtC and by the humus with 0.79 GtC. For the simulation without terrestrial organic matter inputs, the 485 ocean behaves differently and rather uptakes carbon with a slightly negative CO 2 flux. This opposite behaviour between the two simulations with and without terrestrial organic matter fluxes is also observed in central and sSouthern Indonesia and highlights the key role of terrestrial organic matter fluxes in the oceanic outgassing. In central partIndonesia (4.5°N -107.5°E ) we observe an outgassing peak of 8.8 x 10 −9 kg C m −2 s −1 at 14.54 ka and then a positive CO 2 flux for 300 years ( Figure   11b). This outgassing event is primarily the consequence of wood input to the ocean with 1.95 GtC, as well as humus input, 490 the second largest contributor, with 0.69 GtC. In sSouthern partIndonesia (6.5°S -134.5°E) we also have an outgassing with 5.3 x 10 −9 kg C m −2 s −1 to the atmosphere at 14.18 ka due to dominant wood inputs of 2.15 GtC (humus inputs represent 0.67 GtC) (Figure 11c).

Aerobic remineralization
For a higher carbon to nutrients ratio the three locations in Indonesia reproduce similar CO 2 fluxes to the atmosphere than those in the reference simulation. For a lower carbon to nutrients ratio, the CO 2 flux is smaller than the reference simulation 495 with only slightly positive values following the flooding event, with maximum of 0.085 x 10 −9 kg C m −2 s −1 for nNorthern partIndonesia, 0.67 x 10 −9 kg C m −2 s −1 for central partIndonesia and 0.14 x 10 −9 kg C m −2 s −1 for sSouthern partIndonesia ( Figure 11). These sensitivity experiments highlight the fact that even with a very high carbon to nutrients ratio (so larger fraction of remineralized terrestrial carbon in the water column), the outgassing CO 2 flux doesn't increase. However, in the case of the lower carbon to nutrients ratio for terrestrial organic matter, the CO 2 flux to the atmosphere is greatly reduced. But 500 again, these fluxes are only happening at regional hotspots and do not affect the global net CO 2 flux between the ocean and the  Besides the C:N:P ratios, the remineralization rates of the terrestrial organic matter in sea water are not well constrained parameters. The choice of different rates could lead to higher or lower CO 2 flux to the atmosphere. In the deglaciation run and presented sensitivity simulations with higher and lower stoichiometries of terrestrial organic matter, the remineralization rates This simulation uses the same higher stoichiometry ratios as one of the the first sensitivity studies (see Table 3) to get an upper estimate of the potential impact of terrestrial fluxes.
We observe higher CO 2 outgassing in the defined equatorial box over a shorter time period (Figure 11). For the northern part, 510 the CO 2 flux to the atmosphere reaches 20 x 10 −9 kg C m −2 s −1 after the flooding at 14.64 ka and decreases twice as fast as the simulation with high stoichiometry (Figure 11a). Similar behaviour is observed for central and southern part with an outgassing peak after the flooding at 14.54 ka and 14.18 ka of respectively 32 x 10 −9 kg C m −2 s −1 and 27 x 10 −9 kg C m −2 s −1 (Figure 11b,c). This increased CO 2 flux to the atmosphere is primarily a result of wood remineralization since, as for previous simulations, wood dominates the terrestrial organic matter input to the ocean during flooding events at that latitude.

515
Since wood is not buried in the sediment, the amount of material that can be remineralized is the same as in previous simulation, but at faster rate. Part of the outgassing is still due to the remineralization of woody litter and humus before they are buried.

Conclusions
In this study we present for the first time the implementation of terrestrial organic matter fluxes between land and ocean at transiently changing land-sea interface in the global Earth System Model MPI-ESM. This unique setup of MPI-ESM was used 520 to perform a transient deglaciation simulation from 21 to 12 ka accounting for sea level rise induced by meltwater inputs from ice sheets and consequential changes in ocean depth and coast lines. The period between 15-14 ka, which corresponds to Meltwater Pulse 1a, has been highlighted because it is characterized by larger terrestrial organic matter inputs to the ocean than during the first part of the deglaciation. Indeed, a total of 21.2 GtC, mostly arising from Indonesia, goes to the ocean during this millennial event, which represents 34 % of the total amount of terrestrial organic matter entering the ocean over the last 525 deglaciation. This terrestrial carbon is remineralized once in the ocean within a time frame of hundreds of years. The effect of this supplementary carbon brought from land is only observed at regional hotspots with local outgassing to the atmosphere.
The carbon input doesn't seem to be large enough to impact the global behaviour of the ocean, considering that it represents a very small amount in comparison to the global ocean inventory (0.06 %). A sensitivity experiment also emphasizes that the terrestrial organic matter fluxes only have a small effect on the surface alkalinity and dissolved inorganic carbon (around 1-2 530 % increase). However, the local CO 2 fluxes between the ocean and the atmosphere during MWP1a are driven by the terrestrial organic matter inputs. This regional outgassing to the atmosphere observed in Indonesia is explained by wood inputs and is supported by several lines of evidences suggesting that prior to this meltwater pulse event and even as far back as the Last Glacial Maximum, tropical forest was developed in this region favouring the storage of carbon-rich material on land entering the ocean once the land is flooded. As a complement, an additional set of sensitivity experiments show that the magnitude of 535 outgassing during MWP1a is rather insensitive to higher carbon to nutrients ratio of the terrestrial organic matter but rather responds to higher remineralization rates of terrestrial organic matter. Overall, our simulation is a first step towards a fully coupled ESM including carbon and nutrients fluxes at the land-sea continuum that will be applicable for long transient paleosimulations over the last glacial/interglacial cycles. Heinze, C. and Maier-Reimer, E.: The Hamburg oceanic carbon cycle circulation model version "HAMOCC2s" for long time integrations, DKRZ, 1999.