A Persistent Kuroshio in the Glacial East China Sea and Implications for Coral Paleobiogeography

The Kuroshio Current is a major hydrographical feature of the modern East China Sea, but it has been suggested that its flow was diverted to the east of the Ryukyu Arc at the Last Glacial Maximum. Shoaling of the Yonaguni Depression has also been proposed as a cause of Kuroshio Current diversion which, while unlikely to have been significant at the Last Glacial Maximum, may have been an important consideration further back in time. Using an ensemble of high‐resolution ocean simulations with climatic boundary conditions emulating those of the Last Glacial Maximum, we present the first regional state estimates of the glacial East China Sea which are both physically consistent and compatible with sea surface temperature proxy compilations. We find that while the Kuroshio Current transport in the East China Sea is slightly reduced at the Last Glacial Maximum, its path is relatively unchanged, with limited sensitivity to glacioeustatic sea level change, glacial‐interglacial changes in climate, and tectonic shoaling of the Yonaguni Depression. Simulations with the best model‐proxy agreement predict only limited contraction of the reef front at the Last Glacial Maximum, and strong surface currents associated with the glacial Kuroshio may have maintained or even improved long‐distance coral larval dispersal along the Ryukyu Arc, suggesting that conditions may have enabled coral reefs in this region to remain widespread throughout the last glacial. Further field studies investigating whether this is genuinely the case will provide insights into how the coral reef front responds to long‐term environmental change.


Introduction
The Kuroshio Current is the western boundary current of the North Pacific subtropical gyre. Originating as the northward offshoot of the North Equatorial Current, the present-day Kuroshio enters the East China Sea (ECS) through the Yonaguni Depression (YD) east of Taiwan (Figure 1a). It then flows as a topographically steered current along the Okinawa Trough (OT), before leaving the ECS south of Kyushu through the Tokara Strait (e.g., Andres et al., 2015). The vast heat transport associated with the Kuroshio elevates surface temperatures in the ECS and facilitates a significant coral larval flux along the Ryukyu Islands, enabling Japan to host some of the northernmost known coral reefs (Yamano et al., 2011(Yamano et al., , 2012. and present (e.g., Iryu, 2006;Price et al., 2019;Yamano et al., 2011). In particular, the Last Glacial Maximum (LGM) and subsequent recovery therefrom may serve as a useful case study for this purpose. There has therefore been considerable interest in reconstructing coral paleobiogeography and controls such as sea surface temperature (SST) from the glacial ECS. At the LGM, SST was depressed in the ECS (e.g., Kim et al., 2015), and the coral reef front is thought to have retreated along the Ryukyu Arc (Iryu, 2006). Corals of LGM-age have been recovered off Irabu island in the southern Ryukyus (Sasaki et al., 2006), and reef-like structures have been identified in LGM-aged strata in seismic lines from the central Ryukyus (Arai et al., 2012;Matsuda et al., 2011) (Figure 1a).
As the dominant hydrographical feature in the East China Sea, the path and strength of the glacial Kuroshio are inextricably linked to discussions on glacial SST. Based on planktonic foraminifera assemblages, Ujiié and Ujiié (1999) suggested that the Kuroshio was diverted east of the Ryukyu Arc into the Ryukyu Trench due to the formation of a hypothetical Ryukyu-Taiwan land bridge. However, subsequent geological and geochemical support for a displaced glacial Kuroshio is equivocal (e.g., Gallagher et al., 2015), and the existence of such a land bridge at the LGM would require exceptional rates of extension in the Okinawa Trough to explain the subsequent submergence of the YD (see supporting information Text S1). It is however reasonable to expect significant changes in the sill depths (such as the YD) around the OT further back in time (e.g., Koba, 1992;Kao, Roberts, et al. 2006;Lallemand et al., 1999).
Displacement of the Kuroshio at the LGM has been explored from a modeling perspective, testing the effect of glacioeustasy only on the path and strength of the Kuroshio (Kao, Wu, et al. 2006; Lee et al., 2013). Both studies simulated a reduction in the strength of the Kuroshio in the ECS, but while one found that lower sea level caused the Kuroshio to leave the OT via the Kerama Gap instead of the Tokara Strait (Kao, Wu, et al. 2006), the other found little effect of sea level on the Kuroshio path (Lee et al., 2013). These studies were limited to investigating the effects of sea level change, although both acknowledged that other factors such as climate may also be important, which we investigate here. While Lee et al. (2013) also used selected SST  Table 2 (section C corresponds to the PN-line). The contemporary path of the Kuroshio Current is schematically represented by the translucent white arrow. Confirmed LGM-aged reefs at Irabu Island (Sasaki et al., 2006) and suspected LGM-aged reefs off Okinawa and Amami Island (Arai et al., 2012;Matsuda et al., 2011) are plotted as black and white stars, respectively. The suspected LGM-aged reefs have been identified from seismic lines but have not yet been sampled. Also plotted are the core locations of SST proxies used in this study (see Table S2 and references therein), subdivided into proxies based on Mg/Ca (blue squares), alkenones (orange circles), and planktonic foraminifera/radiolarian assemblages (green triangles). (b) Map of the western North Pacific with the model domains used in this study. The outer domain (L2) is forced at the boundaries by PMIP3 climatologies. An extended L2 domain (dashed) is used for pre-industrial (PI) simulations to avoid an open boundary in the Yellow Sea. The L1 domain is nested within the L2 domain with a two-way transfer of information (Debreu et al., 2012). Background is from www.shadedrelief.com. proxies to test the relative effect on SST of forcing the Kuroshio to leave the OT via the Kerama Gap, using realistic climate boundary conditions enables direct comparison with SST proxy compilations. We can thereby quantify model-proxy agreement and determine a best estimate for the state of the glacial ECS. Finally, in addition to impacts on SST, accurate simulations of the glacial Kuroshio will contribute toward understanding coral larval dispersal during the LGM and deglaciation.
Here, we present the first simulations of the glacial Kuroshio under realistic bathymetric and climatic boundary conditions. Forced by multimodel paleoclimate ensembles, enabling direct model-proxy comparisons, we have designed a suite of simulations to answer the following questions: 1. How sensitive is the path and strength of the Kuroshio in the East China Sea to glacial-interglacial changes in climate and glacioeustasy and tectonic shoaling of the Yonaguni Depression? 2. How may changes in SST and current strength have affected coral reef paleobiogeography in the East China Sea at the LGM?

Methods
Experiments were performed using the Coastal and Regional Ocean Community Model (CROCO) V1.0 in a two-way nested configuration (Debreu et al., 2012). This configuration includes 50 vertical s-coordinate layers with increasing resolution at the surface and sea-floor and outer (L2) and inner (L1) domains with resolution 1/16°(∼⃒ 6km) and 1/48°(∼⃒ 2km), respectively ( Figure 1b). These simulations are therefore run at a considerably higher resolution than the 1/8°resolution used in previous studies (Kao, Wu, et al. 2006;Lee et al., 2013). The Yonaguni Depression, where the Kuroshio enters the East China Sea, is narrow. Taiwan and Yonaguni Island, between which the Kuroshio flow is concentrated today, are only separated by about 100km, or around eight grid cells in the previous modeling studies in this region, compared to around 40 grid cells in our highest resolution domain. We therefore used a significantly higher resolution to ensure that the flow through this entrance was adequately resolved, as well as to capture some of the submesoscale features relevant to the dynamics of the Kuroshio (e.g., Kamidaira et al., 2017). Full model parameters can be found in Table S1.
The Paleoclimate Modelling Intercomparison Project Phase 3 (PMIP3) is an attempt to model past climate states through coordinated experiments using global climate models of contributing institutions (Braconnot et al., 2012). Although these models are too coarse to properly resolve regions such as the ECS, the output can be used as a boundary condition for higher resolution models that are able to resolve processes relevant to regional climate (Ludwig et al., 2019). For our experiments, we used monthly climatological pre-industrial (PI) and LGM boundary conditions from four PMIP3 ensemble members, chosen on the basis of resolution and representation of the ECS region: MRI-CGCM3 (Yukimoto et al., 2012), CNRM-CM5 (Voldoire et al., 2013), MIROC-ESM (Watanabe et al., 2011), and GISS E2-R (Schmidt et al., 2014) (Table S1, hereon MRI, CNRM, MIROC, and GISS). References to these ensemble members hereon refer to our simulations driven by the boundary conditions from these models, rather than the original models. All boundary conditions were pre-processed for use in the model using the CROCO_TOOLS V1.0 package, and surface forcing fields were further processed with a Gaussian filter to minimize downscaling artefacts. Rivers are not explicitly included due to insufficient constraints at the LGM. Tides (and paleo-tides) are also not included in these simulations as they are unlikely to be a primary control on the Kuroshio strength (e.g., Ichikawa & Beardsley, 1993). Models were initialized from the interpolated parent PMIP3 outputs and were integrated for 8years, with the final 5 years used to produce the results presentedhere.
Bathymetry is based on the 15 arc-second GEBCO 2019 grid (www.gebco.net). To generate the bathymetry for the LGM runs, we assumed a eustatic sea level fall of 135 m (Lambeck et al., 2014). We do not consider the relative component of glacial isostatic adjustment as this was not large enough in the ECS region to significantly affect any of the major straits (Milne & Mitrovica, 2008), with the possible exception of the Tsushima Strait (Gallagher et al., 2015;Park et al., 2000). The ECS is bounded by an active subduction zone, the Ryukyu Arc, and while significant rates of uplift have been calculated at certain locations (Ota & Omura, 1992), these appear to be geographically restricted to the forearc region (Henry et al., 2014). We therefore consider this an acceptable uncertainty and assume only glacioeustatic sea level change aside from discussion pertaining to the depth of the YD sill.
The PI runs were validated against the World Ocean Atlas 2018 SST climatology (Locarnini et al., 2019) (hereon WOA18) and a reconstructed PI SST climatology using the HadISST1 data set (Rayner et al., 2003). To validate the LGM runs, we principally relied on proxy data compiled by Kim et al. (2015), specifically planktonic foraminiferal Mg/Ca, the alkenone unsaturation index U K ′ 37 , and planktonic foraminiferal assemblages. In their review, Kim et al. (2015) also calculate the seasonalities that these proxies best represent in the OT, which we use here. We direct the reader to this review for a discussion on how these proxies work. We supplemented these data with all other proxy ECS SST estimates that we could find in the literature from the LGM or the earliest deglacial period (Table S2). To directly compare modeled SST to proxy SST, we calculated "simulated" proxy SST by interpolating modeled SST to the core coordinates, using the seasonalities calculated by Kim et al. (2015). Core sites within the L1 domain used the modeled L1 SST, whereas the small number of core sites in the northernmost OT that falls outside of the L1 domain used the modeled L2 SST ( Figure 1a). We also compared modeled and proxy changes in SST between PI and the LGM (ΔSST), with the core top SST used for PI if available or late Holocene SST if not.
The habitable range of coral reefs in the glacial ECS was predicted by extracting the winter minimum 18°C isotherm (see section 4) from the climatological modeled SST. A second estimate of this isotherm (LGM * ) was produced by using the modeled PI-LGM ΔSST instead of the LGM SST field. The modeled PI-LGM temperature change was subtracted from the winter minimum SST from the WOA18 climatology (thus retaining all spatial heterogeneity in ΔSST), producing a new estimate for glacial winter minimum SST. The LGM * winter minimum 18°C isotherm was then extracted from this new SST field.
We have designed three sets of experiments to test the questions posed in section 1, summarized in Table 1. To produce a realistic state estimate for the glacial Kuroshio, we first ran paired PI and LGM simulations for each of the four PMIP3 ensemble members described above to evaluate glacial-interglacial changes in the state of the Kuroshio Current and East China Sea and then quantitatively compared the model output to proxy data. To test the effects of varying YD sill depth, we augmented our LGM simulation with the best model-proxy agreement by synthetically shoaling the maximum depth of the YD from the glacial default to 600, 400, 200, and 0 m and compared the downstream effects on the Kuroshio and ECS to proxy data. Finally, we ran an additional simulation under PI climate and boundary conditions and LGM sea level to decouple these two forcings on the path and strength of the Kuroshio.

Validation
Comparison of the PI SST fields with the WOA 2018 SST climatology reveals differing agreement among the four ensemble members used ( Figure S1 and Table S3), with the model forced by CNRM boundary conditions performing the best. This SST difference is partly due to climate change since the mid-19th century (the time period the PI case simulates). Removing this effect by subtracting the change in SST since 1870 from the 1°gridded HadISST1 data set significantly improves the model-data fit (Table S3). A warm bias remains in the model residual in the southwestern ECS associated with high flow through the Taiwan Strait, and a cold bias exists along the Chinese coast due to anomalous cold water export from the Yellow Sea. Neither the Taiwan Strait nor the Yellow Sea existed at the LGM. These anomalies aside, modeled PI SST is reasonable.
Velocity transects across the PN-line ( Figure S2) demonstrate that all ensemble members reproduce the key features of the Kuroshio in the ECS from observations (Andres et al., 2008) including the depth structure of the Kuroshio itself, the deep countercurrent, and the offshore recirculation. Annual mean net transport across the PN-line varies between 25.8 and 34.3 Sv across the ensemble ( LGM MRI-CGCM3 LGM LGM (CNRM) LGM CNRM-CM5 LGM LGM (MIROC) LGM MIROC-ESM LGM LGM (GISS) LGM GISS E2-R LGM YD shoaling experiments Control (LGM CNRM) LGM CNRM-CM5 LGM 600 m deep YD LGM (YD ≤ 600 m) CNRM-CM5 LGM 400 m deep YD LGM LGM ( , 2017). This overestimate may be driven by a weak recirculation compared to observations (Andres et al., 2008), although it is notable that several other high-resolution models covering the ECS also produce high estimates for transport across the PN-line (Kamidaira et al., 2017;Miyazawa et al., 2009).
Turning to the LGM, Figure 2e compares (geochemical) proxy LGM SST from sites across the ECS with modeled LGM SST. The geochemical proxies (Mg/Ca and alkenones) are generally consistent with modeled LGM SST, whereas those based on foraminifera assemblages perform very poorly, particularly in the winter ( Figure S3a). In their review of SST proxies in the ECS, Kim et al. (2015) reach the same conclusion and argue that "alkenones and Mg/Ca appear to be the most robust choices for paleothermometry in the OT." We have therefore excluded these proxies from comparisons here, but full comparisons are provided in Table S4. Modeled SST is generally within 2°C of the (geochemical) proxy estimates, with the exception of the sites in the northernmost OT. Spread among the models is greatest here due to disagreement regarding the penetration of the Kuroshio into these shallower parts of the OT, and better model-proxy agreement is obtained by models that predict milder temperature reductions in this region. The best performing model is CNRM (which also performed best in the PI case), followed by MIROC, with RMS model-(geochemical) proxy offsets of 1.39°C and 1.48°C, comparable to proxy uncertainties (Kim et al., 2015). Aside from CNRM, most ensemble members generally predict cooler temperatures than the proxies. The poorest performing member is GISS.
An alternative method of validating the LGM models is comparing the modeled change in SST (ΔSST) from PI to the LGM (Figure 3e), which may remove the effect of persistent biases in either the model or proxies. This improves the RMS model-proxy offset for all ensemble members apart from MIROC. The best performing model by this metric is MRI, with its RMS model-proxy offset (1.00°C) almost halved relative to the absolute SST metric (1.93°C). Both CNRM and MIROC predict a smaller temperature drop at the LGM compared to proxies. As with the absolute SST metric, GISS performs the poorest.

The Glacial Kuroshio
The Kuroshio core is often defined as being bounded by the 0.4 m s −1 surface isotach (Andres et al., 2008;Guan, 1980) and is plotted for the all PI-LGM pairs in Figure 4. The LGM Kuroshio core is displaced slightly into the OT relative to PI, a consequence of sea level change as the Kuroshio is topographically steered by the continental slope in the ECS (Andres et al., 2015), which also results in a general narrowing of the Kuroshio core. The departure-point of the Kuroshio from the ECS is also shifted southward, exploiting the deepest point of the Tokara Strait. Annual mean surface velocities averaged across the ensemble are shown for the PI and LGM cases in Figure 5, demonstrating that both the PI and LGM simulations are characterized by strong surface flow throughout the OT. The arc-ward shift in the Kuroshio core can be seen by the reduction in current strength at the LGM paleo-coastline, compensated by an increase in current strength closer to the Ryukyu islands (Figure 5c). Transport through the PN-line is reduced by 13-33% (mean 26%) and is compensated for by a significant increase in transport in the Ryukyu Trench (Table 2), which is also seen in Figure 5c as an increase in surface flow along the Pacific side of the Ryukyu Arc. However, critically, the path of the Kuroshio core is not dramatically altered in any simulation, and the current maintains a significant presence throughout the ECS.
Figures 2a-2d and 3a-3d respectively show the winter minimum SST for LGM experiments and change in winter SST (ΔSST) for PI-LGM pairs. We emphasize winter SST due to its role in coral habitability (see section 4). Here, differences between ensemble members become more apparent. Winter minimum SST in the OT is warmest in the CNRM and MIROC cases, falling by only 1-3°C at the LGM relative to PI. The northernmost OT is an exception: The slight southward shift in the departure point of the Kuroshio understates the relatively significant reduction in the influence of the warm Kuroshio waters Note. All transport values refer to the net annual mean transport across that section. Please see Figure 1a for section locations and Table 1 for experiment details 10.1029/2020PA003902

Paleoceanography and Paleoclimatology
VOGT-VINCENT AND MITARAI LGM-aged reefs are plotted as black and white stars, respectively, as in Figure 1a. (e) Offset between modeled and proxy LGM SST for the four LGM experiments (geochemical proxies only: blue squares = Mg/Ca, orange circles = alkenones). The width of the gray envelope represents data-point density. Perfect model-proxy agreement for LGM SST is represented by the horizontal line. Points plotting above the line represent warmer modeled SST than proxy SST. Points plotting below the line represent cooler modeled SST than proxy SST. Numbers in gray represent the RMS model-proxy offset.
in the shallower parts of the OT southwest of Kyushu. Consequently, SST here falls in excess of 3°C. Compared to CNRM and MIROC, MRI predicts somewhat greater cooling, with a temperature drop of 3-4°C throughout most of the OT. GISS is an extreme case, simulating a temperature drop of more than 4°C in the OT.

Sensitivity to the Yonaguni Depression Sill Depth
While there is no good evidence for the exceptional subsidence rates demanded by a land bridge at the LGM as has been previously suggested (Ujiié & Ujiié, 1999), combinations of tectonics and glacioeustasy could conceivably significantly shoal the depth of the YD further back in the geological record. However, our Yonaguni Depression shoaling experiments demonstrate that the Kuroshio in the ECS is remarkably resilient in the face of this shoaling. Even with a sill depth of just 200 m, while weakened, the Kuroshio core still remains entirely within the OT (Figure 6), and changes in SST rarely exceed 0.5°C relative to the unshoaled sill (Figure 7). When the sill is eliminated entirely, the Kuroshio instead enters the OT east of Ishigaki Island and through the Kerama Gap and transport through the PN-line is only reduced by 36% relative to the control (Table 2).
In the reference LGM CNRM simulation, flow through the YD is concentrated in the upper 400 m and to the west of Yonaguni Island, as is the case today (Andres et al., 2015). Indeed, flow through the YD only appears to respond significantly to tectonic shoaling once it becomes shallower than this depth ( Figure S4). When the YD is shoaled to 400 m, transport through the YD is still strong (Table 2), but the depth structure becomes much more uniform, and the peak velocity migrates from the surface to the subsurface. When the YD is shoaled to 200 m, flow through the YD is unsurprisingly significantly suppressed (approximately halving relative to the LGM CNRM reference) and the flow becomes strongly depth-intensified, resulting in a significant localized reduction in mean surface velocity. The relative importance of flow east of Yonaguni Island also increases, although the bulk of transport still occurs in the west. In all cases, however, the structure of the Kuroshio by the PN-line is not significantly affected aside from a general reduction in transport. As a result, with the exception of the southernmost OT, the effect of a Taiwan-Ryukyu land bridge on SST is generally <1°C, small compared to glacial-interglacial changes.

Decoupling Sea Level Change and Climate
While previous modeling studies have focused on the effects of changing sea level on the Kuroshio, it is acknowledged that this may not be the only relevant parameter (Kao, Wu, et al., 2006;Lee et al., 2013). To decouple the effects of sea level change and climate, we have run a simulation using PI CNRM surface and lateral boundary conditions and LGM bathymetry. The effects on transport through the sections (Table 2) and the Kuroshio path in the OT are very similar to the full LGM CNRM simulation, validating the assumption that sea level is the main driver of glacial-interglacial variability of the Kuroshio strength in the ECS (Kao, Wu, et al. 2006;Lee et al., 2013). The location of the Kuroshio core within the ECS is essentially unaltered in this simulation compared to the full LGM simulation although the core axis is more sinuous and displaced eastward southeast of Taiwan ( Figure S5). For SST, climate generally dominates the signal, but bathymetry is locally important (Figure 8). The slight southward shift in the departure point of the Kuroshio and reduction in influence of the Kuroshio in the northernmost OT results in a significant reduction in SST off the southwest coast of Kyushu, similar in magnitude to the effect of climate (Figures 8a and 8c). The slight displacement of the Kuroshio core into the OT also results in warming along the Ryukyu Arc, partially compensating for the climate-induced cooling. Elsewhere in the OT and in the Ryukyu Trench, the effects of bathymetry on SST are secondary to climate (Figures 8b and 8c).

Discussion and Conclusions
We have produced the first estimates for the state of the glacial ECS which are both demonstrably physically consistent and compatible with SST proxy compilations. It is important to be cautious of over-interpreting these results: There is significant uncertainty in the PMIP3 ensemble (Harrison et al., 2015(Harrison et al., , 2016, and local forcings such as changing freshwater fluxes (Ijiri et al., 2005;Xu & Oda, 1999) are poorly constrained and not explicitly incorporated. Tides were also not included in this study, although sea level fall at the LGM resulted in significant changes in tidal dissipation and amplitudes relative to the present (Wilmes & Green, 2014). As a component of the wind-driven gyre circulation, tides are unlikely to be of first-order importance for the strength of the Kuroshio (e.g., Ichikawa & Beardsley, 1993), so the exclusion of tides from our simulations is unlikely to significantly alter our results. However, this should be investigated in future work, particularly pertaining to paleo-dispersal pathways as discussed below. Furthermore, due to the very significant downscaling in these models versus the resolution of the boundary conditions used, it is inevitable that neither the PI nor LGM ensembles simulated here will be a perfect reflection of reality. However, the good model-proxy agreement for both LGM SST and PI-LGM ΔSST, directly comparable to proxy uncertainty for the best ensemble members, indicates that our simulations are reasonable. Significantly, the Kuroshio persists as a major feature in the ECS in all cases, even in the face of significant shoaling or elimination of the YD. We therefore argue that this result is robust.
Which of the four ensemble members represents the best estimate for the state of the glacial ECS? The best model-proxy agreement for absolute LGM SST is found in CNRM and MIROC (Figure 2e), both of which simulated the smallest reduction of SST in the OT (1-3°C). In both of these ensemble members, the smallest temperature changes are on the island-facing side of the OT. This is likely due to the particularly strong increase in current strength on the island-facing side of the OT that both members predict for the LGM ( Figure S6). However, the ΔSST metric, which may help remove any consistent biases in both the model

Paleoceanography and Paleoclimatology
VOGT-VINCENT AND MITARAI and proxies, indicates that MRI, which simulates more cooling in the OT (3-4°C), has the best model-proxy agreement (Figure 3e). One caveat for the ΔSST metric is that the proxy PI SST values are actually a mixture of core top and late Holocene SST. For references that reported both core top and late Holocene SST, some late Holocene values were slightly higher than the core top value. The proxy ΔSST value may therefore be a slight overestimate, which could improve the agreement of CNRM and MIROC. We would therefore conclude that MRI, CNRM, and MIROC all provide useful information about the state of the glacial ECS.
Our results more generally suggest that the impacts on the Kuroshio and surface conditions in the ECS of changing bathymetry, either from eustatic sea level change or tectonic shoaling of the YD, are limited. Comparison of our realistic state estimates to simulations of the effect of eustatic sea level change only ( Figure 8) suggests that glacial-interglacial variability in SST in the ECS is dominated by the climatic signal. The exception is the northernmost OT, where Kuroshio penetration is significantly suppressed by the sea level fall due to the southward shift in the Kuroshio core as it leaves the ECS through the Tokara Strait. The quality of SST predictions for the ECS is therefore strongly dependent on the quality of the climatic boundary conditions, emphasizing the importance of reducing uncertainty in future PMIP contributions to improve the ability of regional paleoceanographic models to predict past SST. The effect of tectonic shoaling of the YD on the Kuroshio also appears to be relatively minor until the sill depth (specifically west of Yonaguni Island) shoals above 400 m. Even when the YD is closed entirely (the hypothetical Taiwan-Ryukyu land bridge), as long as the Kerama Gap remains open, the Kuroshio remains in the central and northern OT. Consequently, the ability of SST proxies to descern the existence of significant shoaling of the YD is limited. Indeed, due to the absence of good LGM SST estimates from the southernmost OT, SST proxies by themselves are currently unable to clearly distinguish between the presence or absence of a hypothetical Taiwan-Ryukyu land bridge ( Figure S3b). These predictions may help guide future proxy studies attempting to reconstruct environmental change in the ECS in response to long-term tectonic changes (e.g., Kimura, 2000). It would also be interesting to explore the application of new proxies for past Kuroshio strength such as the geostrophic method using oxygen isotopes (Lynch-Stieglitz et al., 1999) as reconstructions based on provenance analyses and sortable silt have resulted in divergent interpretations (e.g., Diekmann et al., 2008;Dou et al., 2010;Zheng et al., 2016).
Today, the winter minimum 18°C isotherm approximates the northern limit of coral reefs in Japan (Ministry of the Environment & Japanese Coral Reef Society, 2004), and while this is not the strict limit for coral habitability (e.g., Yamano et al., 2012), it is generally considered to demarcate the limit for reef formation in the region (Veron & Minchin, 1992). When considering changes in the coral habitable range in past climates, it is therefore reasonable to consider the migration of this isotherm as a first approximation. This isotherm is plotted in Figure 2a-2d for the four LGM state estimates we have produced. Figure 3a-3d also show a second isotherm (LGM * ), which assumes that the modeled PI-LGM ΔSST is more accurate than the modeled absolute SST, subtracting ΔSST from the WOA18 climatology. Particularly in the MRI case, this is likely a more accurate isotherm due to the significantly better model-proxy agreement for ΔSST relative to absolute SST, with the smallest model-proxy RMS offset of any comparison in this study (1.00°C).
CNRM and MIROC predict that the winter minimum 18°C isotherm is displaced only moderately across the Tokara Strait, agreeing well with the present distribution of confirmed and suspected LGM-aged coral reefs (Figure 2b-2c). The MRI LGM * isotherm also places this isotherm in the northern OT (Figure 3a). Conversely, this isotherm all but leaves the ECS for GISS and the regular MRI case. This is not conclusively contradicted by the coral data as firstly, modern coral reefs are exceptionally found north of the winter minimum 18°C isotherm (Yamano et al., 2012) and secondly, only the southernmost LGM-aged coral reef (which is relatively close to the GISS and regular MRI winter minimum 18°C isotherms) has been sampled and therefore proven to be as such. However, since all isotherms with the best model-proxy agreement (CNRM, MIROC, and MRI LGM * ) remain in the northern OT, a limited retreat of the winter minimum 18°C isotherm is more compatible with the proxy evidence. This adds significant credibility to the structures identified as reefs from seismic lines by Matsuda et al. (2011) and Arai et al. (2012) and generates the prediction that LGM-aged reefs could conceivably be ubiquitous across much of the Ryukyu Arc. If the site identified by Matsuda et al. (2011) represented the northernmost extent of reef formation in the glacial ECS, this would further indicate that the reef front tracks the winter minimum 18°C isotherm in response to long-term climate change. Alternatively, given that this isotherm is not a strict limit to reef formation (Yamano et al., 2012), it is conceivable that LGM-aged reefs could be discovered further north still. However, for now the extent of reef front retreat at the LGM remains poorly constrained and sampling is required from the sites identified from seismic lines to confirm their age and lithology.
While temperature is the dominant control on the coral habitable range (Veron & Minchin, 1992), other relevant parameters include salinity, turbidity, nutrients, and carbonate saturation state, which are not discussed here but have been simplistically included in modeling studies (Kleypas, 1997). Understanding these controls is important as ocean acidification in particular may damage the potential of temperate oceans to act as coral refugia (Van Hooidonk et al., 2014). Unfortunately, these parameters are at present very challenging to validate and realistically model for past climates. Furthermore, caution must be taken when extrapolating coral reef sensitivity from the late Pleistocene as the warmer climate of the Anthropocene may make relevant new limitations on reef expansion such as the latitudinal attenuation of sunlight (Muir et al., 2015).
Aside from the effects on SST, another important consideration for the effects of interglacial-glacial Kuroshio changes on coral reefs in the ECS (and more generally, the regional paleobiogeography of organisms which disperse through ocean currents) is the impact on propagule dispersal. While not explored in detail here, the model resolves significant changes in the wakes of islands in the path of the Kuroshio (particularly the Yaeyama and Tokara islands) in the LGM simulations. These changes in wake behavior are likely related to changes in the flow past ocean islands caused by the shifts in the Kuroshio core axis (e.g., Heywood et al., 1996) and could affect the retention or dispersal of coral larvae (e.g., Largier et al., 2004;Wolanski & Kingsford, 2014).
Likely more important, however, is the role of the Kuroshio in the long-distance dispersal of coral larvae along the Ryukyu Arc (Uchiyama et al., 2018;Yasuda et al., 2009). While the transport of the Kuroshio was likely reduced at the LGM, this is primarily related to spatial confinement of the current within the OT due to the exposure of the continental shelf, resulting in a narrowing of the Kuroshio core. In our simulations, mean surface current strength within the Kuroshio core was not significantly different at the LGM (0.62 m s −1 ) than for the PI case (0.61 m s −1 ). Furthermore, a present bottleneck for long-distance Kuroshio-mediated dispersal is the entrainment and release of propagules into and from the Kuroshio due to the spatial separation between the Kuroshio core and Ryukyu islands. In all simulations, there was a significant increase in surface velocities associated with the Kuroshio on the Ryukyu-facing side of the Kuroshio axis ( Figures 5 and S6), potentially indicating that long-distance propagule dispersal through the Kuroshio may in fact have been more effective at the LGM than today. The continued presence of strong, unidirectional surface currents along the Ryukyu Arc associated with the Kuroshio at the LGM may therefore have contributed to the resilience of reefs living at their habitable limit during this time of environmental stress and tentatively may have supported the existence of yet-undiscovered glacial-aged coral communities living beyond the 18°C isotherm, as is the case today (Yamano et al., 2012). The persistence of the Kuroshio in the ECS may also have facilitated the recolonization of the northernmost reef sites occupied today during deglaciation.