Skip to main content

Advertisement

Log in

Global trends of ocean CO2 sink and ocean acidification: an observation-based reconstruction of surface ocean inorganic carbon variables

  • Original Article
  • Published:
Journal of Oceanography Aims and scope Submit manuscript

Abstract

Ocean acidification is likely to impact marine ecosystems and human societies adversely and is a carbon cycle issue of great concern. Projecting the degree of ocean acidification and the carbon-climate feedback will require understanding the current status, variability, and trends of ocean inorganic carbon system variables and the ocean carbon sink. With this goal in mind, we reconstructed total alkalinity (TA), dissolved inorganic carbon (DIC), CO2 partial pressure (pCO2sea), sea–air CO2 flux, pH, and aragonite saturation state (Ωarg) for the global ocean based on measurements of pCO2sea and TA. We used a multiple linear regression approach to derive relationships to explain TA and DIC and obtained monthly 1° × 1° gridded values of TA and DIC for the period 1993–2018. These data were converted to pCO2sea, pH, and Ωarg, and monthly sea-air CO2 fluxes were obtained in combination with atmospheric CO2. Mean annual sea–air CO2 flux and its rate of change were estimated to be − 2.0 ± 0.5 PgC year−1 and − 0.3 (PgC year−1) decade−1, respectively. Our analysis revealed that oceanic CO2 uptake decreased during the 1990s and has been increasing since 2000. Our estimate of the globally averaged rate of pH change, − 0.0181 ± 0.0001 decade−1, was consistent with that expected from the trend of atmospheric CO2 growth. However, rates of decline of pH were relatively slow in the Southern Ocean (− 0.0165 ± 0.0001·decade−1) and in the western equatorial Pacific (− 0.0148 ± 0.0002·decade−1). Our estimate of the globally averaged rate of pH change can be used to verify Indicator 14.3.1 of Sustainable Development Goals.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Institutional subscriptions

Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6
Fig. 7
Fig. 8

Similar content being viewed by others

References

Download references

Acknowledgements

Gridded data of inorganic carbon variables and ocean CO2 uptake calculated in this study are available at the JMA web site: https://www.data.jma.go.jp/gmd/kaiyou/english/co2_flux/co2_flux_data_en.html. The Surface Ocean CO2 Atlas (SOCAT) is an international effort, endorsed by the International Ocean Carbon Coordination Project (IOCCP), the Surface Ocean Lower Atmosphere Study (SOLAS) and the Integrated Marine Biosphere Research (IMBeR) program, to deliver a uniformly quality-controlled surface ocean CO2 database. SOCAT data sets are available at https://www.socat.info. We thank the many researchers and funding agencies responsible for the collection of data and quality control for their contributions to SOCAT. Data of GLODAPv2 are also available at https://www.nodc.noaa.gov/ocads/oceans/GLODAPv2/.We also thank the GLODAP v2 team and data originators. GlobColour data (https://globcolour.info) used in this study have been developed, validated, and distributed by ACRI-ST, France. The altimeter products were produced by Ssalto/Duacs and distributed by Aviso+, with support from Cnes (https://www.aviso.altimetry.fr). This study was supported by MRI's research fund C4 for the studies of ocean biogeochemistry and ocean acidification.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Yosuke Iida.

Appendices

Appendix 1: Derivation of equations to estimate surface ocean TA

Here we describe processes used to derive surface ocean TA from SSS and SSDH by using MLR in accord with the method of Takatani et al. (2014). As described in Takatani et al. (2014), there is a distinct negative correlation between nTA and SSDH, with the exception of some scattered data in a narrow SSDH range that reflect SSS variations (Fig. 

Fig. 9
figure 9

Relationships of nTA versus a SSS and b SSDH in the global ocean. Plots are colored by zone. Abbreviations of zone names are defined in the text

9). We therefore used SSS and SSDH as explanatory variables in the MLR.

$$\mathrm{nTA}=f\left(\mathrm{SSS},\mathrm{SSDH}\right).$$
(8)

We used definitions of zones that were similar but slightly modified from the definitions used by Takatani et al. (2014) for the Pacific Ocean as follows. (1) The high SSDH area (SSDH > 0.8 m) was classified as zone 1WP. This zone represented the northern and southern subtropical gyres and the western equatorial warm pool. The nTAs show a nearly constant value of ~ 2300 μmol·kg−1 because few processes that change nTA occur in this area. (2) The area characterized by SSDHs of 0.2–0.8 m south of the salinity maximum around 25° N was assigned to zone 2EP, whereas the method of Takatani et al (2014) divided this area meridionally along 120° W. This zone represents upwelling areas of the eastern equatorials, off Peru and off Costa Rica. The higher nTAs in 2EP than in 1WP reflect the high nTA of upwelled seawater. (3) The area characterized by SSDHs of 0.2–0.8 m north of the salinity maximum represents the subtropical-subarctic transition and was named 3NP. This zone is characterized by seawater with a wide range of SSSs, from ~ 35, which is representative of subtropical seawater, to around 32 − 33, which is representative of subarctic seawater. The regression explanatory variables therefore included SSS. (4) The area characterized by SSDHs less than 0.2 m was assigned to the North Pacific subarctic gyre (4WSAG) and was characterized by a high and nearly constant value of nTA.

An approach similar to that used in the Pacific was applied to the Atlantic. First, an equatorial area characterized by SSS < 35.5 due to the influence of large rivers such as the Amazon and Congo was designated zone 5EA. The area with SSDH > 0.05 m, which is representative of the tropical and subtropical Atlantic, was classified as 6TA. The use of SSS for zone 5EA could better account for the non-zero nTA source of the large rivers (Lefèvre et al. 2010). The nTA in 6TA was nearly constant and similar to the nTA of zone 1WP in the Pacific. The area characterized by − 0.2 > SSDH > 0.2 m represented the subtropical-subarctic transition, and a northern area characterized by SSDHs less than 0.2 m represented the North Atlantic subarctic gyre; these two areas were classified as 7NEA and 8NA, respectively. In these four zones, SSS and/or SSDH were used for the regressions.

In the Indian Ocean, the Bay of Bengal (9BOB), where there is a large input of freshwater from the Ganges–Brahmaputra River, was distinguished from the rest of the Indian Ocean (10IND) by SSS < 34. The SSS in 9BOB was used to correct the nTA in 9BOB for the influence of the riverine input of nTA (Bates et al. 2006) in the same way that SSS was used to correct the nTA in 5EA in the Atlantic. The nTA in 10IND was nearly constant, despite the inclusion of observations during the southwestern monsoon season, which could have affected the carbonate system in this zone (Sreeush et al. 2019).

The Southern Ocean was distinguished by SSDH < 0.05 m from the Atlantic, by SSDH < 0.65 m from the Indian, and by SSDH < − 0.8 m from the Pacific. Midorikawa et al. (2012) used the 5 °C SST isotherm as a proxy for the subarctic front to estimate TA. In this study, the northern and southern parts of the Southern Ocean were distinguished by SSDH = − 0.1 m as a proxy for the frontal area, and the northern was further divided into three ocean sectors (11SOA, 12SOI, 13SOP, and 14SSO). In these zones, both SSDH and SSS ranged widely and were used as MLR variables.

Table 4 summarizes the definitions of zones and subzones, and the derived regression coefficients for each subzone are listed in Table 5.

Table 4 Zone and subzone definitions used to estimate nTA and nDIC. Abbreviations of zones and subzones can be found in the text
Table 5 Regression coefficients for estimating nTA. Abbreviations of zones are defined in the text

Appendix 2: Derivation of equations to estimate surface ocean DIC

Here, we describe the derivation of the MLR equations used to estimate surface ocean DIC from SST, SSS, SSDH, Chl, and MLD. We particularly emphasize the procedure for defining subzones.

2.1 Data procedure

We used SOCAT V2019 as a surface carbon dataset. Before conducting the MLR analysis, sub-hour data were binned to hourly data to eliminate possible bias from sampling frequency. Some records show unprecedented values of salinity during part or all of each cruise, probably because of equipment malfunction, and we excluded those records from the analysis. After that, fCO2 data were converted to DIC in combination with the observed SST, SSS, the TA estimated from the observed SSS and satellite SSDH data, and climatological nutrients from the World Ocean Atlas 2018 (Garcia et al. 2018). To ignore DIC change resulted from simple fresh water flux, nDIC was used for MLR.

$$\mathrm{nDIC}=f\left(\mathrm{time},\mathrm{SST},\mathrm{SSS},\mathrm{SSDH},\mathrm{MLD},\mathrm{Chl}\right),$$
(9)

where the function f took a more general form than a simple MLR (vide infra). If a large number of explanatory variables are included in a MLR, the coefficients associated with some of the variables will be statistically insignificant. However, this statistical insignificance might not adversely affect the results if the MLR was applied within a limited area where the ranges of the variables were small.

In general, nDIC values were negatively correlated with SST, SSS, and SSDH over most of the ocean (Figs.

Fig. 10
figure 10

Relationships of nDIC versus SST (left), SSS (center), or SSDH (right) in representative subzones of the Pacific Ocean. Coloring indicates the number of data in each grid. Abbreviations of subzone names are defined in the text

10,

Fig. 11
figure 11

The same as Fig. 10, but for the Atlantic Ocean

11,

Fig. 12
figure 12

The same as Fig. 10, but for the Indian Ocean

12 and

Fig. 13
figure 13

The same as Fig. 10, but for the Southern Ocean

13). Based on the relationships in Figs. 10, 11, 12 and 13, we hypothesized that globally nDIC was a quadratic function of SST, SSS and SSDH. In a few subzones, the powers of the explanatory variables were tuned to avoid unlikely extrapolations. We assumed that the secular trends of nDIC were linear. In addition to SST, SSS, and SSDH, Chl and MLD were considered as additional variables based on the a priori hypothesis that Chl would be negatively correlated with nDIC because primary production would be expected to reduce nDIC, particularly from spring to autumn in subpolar subzones. Deepening of the mixed layer during the winter in subtropical and subpolar subzones would increase nDIC because of entrainment from lower seawaters. In some regions, the Chl concentrations were indeed negatively correlated with the nDIC residuals that remained after elimination of the variability associated with SST, SSS, and SSDH (Figs. 

Fig. 14
figure 14

Relationships of nDIC residuals versus SSS (left), Chl (center), and MLD (right) in representative subzones of the Pacific Ocean. Coloring indicates number of data in each grid. Residuals for SSS were calculated via MLR using SST alone as the explanatory variable; those for Chl via MLR using SST, SSS, and SSDH as explanatory variables; and those for MLD via MLR using SST, SSS, SSDH, and Chl as explanatory variables. The x-axes for Chl and MLD are logarithmic. Abbreviations of subzone names are defined in the text

14B-b,

Fig. 15
figure 15

The same as Fig. 14, but for the Atlantic Ocean

15B-b,

Fig. 16
figure 16

The same as Fig. 14, but for the Indian Ocean

17C-b, D-b). The MLDs also showed distinct relationships with the residuals (Figs. 15B-c;

Fig. 17
figure 17

The same as Fig. 14, but for the Southern Ocean

16B-c). We included Chl and MLD in the equations in the form logarithm (c1*X + c2*logX) following Chierich et al. (2012). The following sections explain the classification of subzones in the oceans. The subzones are defined in Table 4, and the derived regression coefficients are listed in Tables

Table 6 Regression coefficients for estimating nDIC. Abbreviations of subzones are defined in the text

6 and

Table 7 Regression coefficients for estimating nDIC in the subtropical Pacific, Atlantic and Indian Oceans. Abbreviations of subzones are written in the text

7.

2.2 Pacific ocean

The Pacific Ocean was divided into four zones for purposes of TA estimation as described in Sect. A1. We further divided it into 16 time-varying subzones, with borders defined on the basis of the distributions of SST, SSS, and SSDH, partly in accord with the methods of Sugimoto et al. (2012) and Iida et al. (2015).

2.2.1 Western equatorial-subtropical subzones (TA zone 1WP)

TA zone 1WP (SSDH > 0.8 m) covered a wide range in the western to central Pacific. The nTA within 1WP was almost constant ~ 2300 μmol·kg−1. At subtropical latitudes in zone 1WP, the carbonate system closely follows a seasonal cycle (e.g., Chierici et al. 2006), whereas non-seasonal variability associated with ENSO characterizes the carbonate system near the equator (e.g., Feely et al. 2006). Based on this difference, zone 1WP was divided into four subzones: two north and south subtropical subzones (1WP-1NWP and 1WP-2SWP, respectively) and two equatorial subzones (1WP-3WEQ and 1WP-4CEQ). At the equator, high salinity, low temperature, and high DIC water is brought to the surface from below by upwelling. There is a salinity gradient along the equator from high salinity in the east to low salinity in the west (Ishii et al. 2009). This part of 1WP was therefore divided into two equatorial subzones, 1WP-3WEQ, the western warm pool with SSS < 34.7 and SST > 27.5, and 1WP-4CEQ, the equatorial divergence area. The boundaries between the subtropical and equatorial regions of 1WP were defined by SSDH minima as proxies for the North and South Equatorial Currents.

Whereas the subtropics were previously thought to be less productive because of a lack of nutrients, it has recently been pointed out that nitrogen fixation and dissolved organic nutrients are important sources of nitrogen for primary production in the ocean (e.g., Hashihama et al. 2013). In addition to biological processes, winter vertical mixing and seasonal expansion and contraction of the subtropical gyres affect nDIC seasonality (e.g., Ayers and Lozier 2012; Fassbender et al. 2017). Consequently, nDICs in summer and winter are different under the same SST and SSS conditions, especially at relatively high latitudes in the subtropics. We took this seasonality into account by formulating monthly equations for the subtropical North Pacific to avoid amplification of the standard errors of the MLR by seasonal effects. We used a similar procedure for the South Pacific, but we used seasonal rather than monthly equations because we had insufficient data for the latter.

2.2.2 Eastern equatorial subzones (TA zone 2EP)

TA zone 2EP, the area with SSDHs of 0.2–0.8 m, was divided into northern and southern parts by the SSDH maximum to the northern rim of the equatorial upwelling zone, and we further divided the two parts into two and three subzones, respectively. The southern three subzones were (1) the central-to-eastern equatorial divergence zone (2EP-3EEQ) defined as SSS > 34 and characterized by high DIC, high salinity, and low temperature (e.g., Feely et al. 1999); (2) the coastal upwelling area off Peru, south of 15°S (2EP-1PERU), which has oceanographic characteristics similar to equatorial (Friederich et al. 2008); and (3) the permanently low-salinity (SSS < 34) area off Panama (2EP-2EPFP) referred to as the eastern Pacific fresh pool (Alory et al. 2012), where both the pCO2 and salinity are low (Brown et al. 2015). The highest nDIC residuals after elimination of SST effects in the 2EP-3EEQ subzone were observed at an SSS of about 34.9 (Fig. 14D-a), which corresponds to a subsurface salinity maximum. The northern two subzones were off Costa Rica (2EP-4COS) to the area east of 105° W and the SSDH less than 0.7 m, which shows the typically high pCO2 related to high wind speed (Brown et al. 2015), and the remainder of the 2EP zone, the northeastern equatorial subzone (2EP-5NEEQ).

2.2.3 Subpolar subzones (TA zones 3NP and 4WSAG)

TA zone 4WSAG, the western subarctic gyre, is characterized by a SSDH less than 0.2 m. The nDIC in the gyre are homogeneous relative to the transition area as the case of nTA, and it was not divided into any subzones during the winter (4WSAG-1w). During the summer, however, it was divided into two subzones (4WSAG-2LTs and 4WSAG-3HTs) at an SST of 7.5 °C, which corresponds to a slight inflection point of the nDIC-SST relationship.

TA zone 3NP, the subtropical-subpolar transition zone, is characterized by steep gradients of SSS and TA as well as SSDH. Seawater in this zone can be considered to be hypothetical mixtures of subtropical and subarctic seawaters (Kakehi et al. 2017). We first defined the southern area of the 3NP zone as a nutrient-depleted, subtropical-like area with SST > 16 °C, and the northern area as a low-SST region around the western subarctic gyre (3NP-6NSBA) with a SST < 5 °C. The subtropical-like area was divided into two subzones by the SSDH maximum along 120° W, an eastern boundary subzone off California (3NP-2CAL), which is affected by CO2 rich subsurface seawater (Alin et al. 2012), and a northern subtropical Pacific subzone (3NP-1NSBT). The intermediate SST (5 < SST < 16 °C) area was divided primarily on the basis of differences between nDIC–salinity relationships into three subzones, a southern subzone with SSDH > 0.35 m (3NP-3STR), a subzone off Sanriku with SSDH < 0.35 m and SSS > 33.3 (3NP-4SRK), and an eastern subzone with SSDH > 0.35 m and SSS < 33.3 (3NP-5ESBA). Overall, in the transition-subarctic subzones, the nDIC residuals after elimination of SST variability were strongly correlated with SSS (Fig. 14C-a), and the nDIC residuals after elimination of SST/SSS/SSDH variability were correlated with Chl concentrations (Fig. 14B-b). Because Chl data are unavailable during winter in high latitudes, these subzones were treated separately during the winter (November–February) and summer (March–October) seasons.

2.3 Atlantic Ocean

2.3.1 Equatorial and subtropical subzones (TA zones 5EA 6TA)

In the equatorial Atlantic, large rivers such as the Congo and Amazon have a large effect on both the salinity and carbonate system (Lefevre et al. 2010; Ibánhez et al. 2015). TA zone 5EA, the area characterized by SSS < 35.5, was defined as a single subzone in the equatorial Atlantic (5EA-1LS). After elimination of SST variability, nDIC residuals in this subzone were strongly correlated with SSS (Fig. 15E-a). TA zone 6TA was divided into eight subzones in a manner similar to the division of Pacific. The western and eastern Atlantic were distinguished by SSDH = 0.4 m. The western Atlantic corresponded to TA zone 1WP and was divided into a subtropical North Atlantic (6TA-1NSBT), South Atlantic (6TA-2SSBT) and western equatorial (6TA-3WEQ) subzones on the basis of SSDH minima in a manner similar to the division of subzones in TA zone 1WP. To disentangle the complexity of the Atlantic eastern boundary systems, we divided the eastern Atlantic into six subzones based on SSDH maxima. Coastal upwelling controls the surface CO2 in the Benguela Current system (Santana-Casiano et al. 2009) and the Guinea dome (Lefevre et al. 1998), and seasonal upwelling lowers SST and increases surface DIC in the eastern equatorial Atlantic (Parard et al. 2010). We distinguished the Angola-Benguela subzone (6TA-4ABR) based on the southeastern Atlantic SSDH maximum along 8° S. The subzone off Guinea (6TA-7GUI) was distinguished on the basis of the SSDH maximum around 15–20° N. The eastern equatorial upwelling subzone (6TA-5EEQ) was located between subzones 6TA-4ABR and 6TA-7GUI. The remainder of the northern area of TA zone 6TA was divided into three subzones, a Gulf Stream subzone (6TA-6GST) west of 60° W; a northeastern subtropical subzone (6TA-8NESBT) defined by SST > 18 °C, a temperature that corresponds to a threshold for nutrient depletion; and a northern subtropics subzone (6TA-9NSBT) defined by SST < 18 °C.

2.3.2 Subarctic subzones (TA zones 7NEA and 8NA)

In the northern North Atlantic, subtropical seawaters are brought to high latitudes by the North Atlantic Current, and seawater characterized by high SST and high SSS is widely distributed in this area (e.g., Rossby 1996). However, low SST and low SSS areas that are influenced by Arctic seawater are distributed around Greenland, the Labrador Sea, and off Newfoundland (e.g., Mertens et al. 2014). These differences of SST and SSS accentuate the contrast of carbonate biogeochemistry between the eastern and western parts of the northern North Atlantic (e.g., Körtzinger et al. 2008; Schuster et al. 2009). The lowering of SST due to deep convection in winter results in a positive correlation between MLD and DIC (e.g., Olsen et al. 2008).

Nutrient-depleted conditions were associated with SST > 18 °C in the subtropics-like subzone (7NEA-1SSBA). These conditions (SST > 18 °C) were found in both TA zones 7NEA and 8NA. The remainder of TA zone 7NEA was defined as a single subzone (7NEA-2SESBA). The remainder of TA zone 8NA was divided into a western subzone with SSS < 34.3 (8NA-1GLB), which mainly included areas around Greenland, the Labrador Sea, and off Newfoundland, and a North Atlantic Current area with SSS > 34.3. The latter was further divided into a southern subzone (8NA-2ESBA) and a northern subzone (8NA-3NSBA) by SSTs greater than or less than 7.5 °C, respectively. Those subzones were distinguished between winter (Nov.–Feb.) and summer (Mar.–Oct.) based on the availability of Chl data. The fact that residuals of nDIC after elimination of effects due to SST, SSS and SSDH were well correlated with MLD as well as Chl (Fig. 15B-c, B-b) indicated that MLD reduced the uncertainty of TA in these subzones.

2.4 Indian Ocean

In the northern Indian Ocean, the monsoons influence the carbonate system in the Arabian Sea, and the Ganges–Brahmaputra River controls the oceanographic system in the Bay of Bengal (e.g., Sabine et al. 2000; Sarma 2003; Bates et al. 2006). The northern Indian Ocean is thought to be a small source of atmospheric CO2, whereas the southern Indian Ocean is a strong sink and absorbs CO2 during the austral winter (Sarma et al. 2013). TA zone 9BOB, the Bay of Bengal, was defined as a single subzone (9BOB-1BOB) in which nDIC was strongly controlled by SSS variability (Fig. 16E-a). In the Arabian Sea and off Somalia, the southwestern monsoon causes coastal upwelling and brings cold seawater to the surface. The upwelling results in surface waters characterized by low temperatures and an extremely high content of CO2 (e.g., Körzinger et al. 1997). In winter, the northeastern monsoon induces convective mixing (Sarma 2003) that results in a high correlation between SST and nDIC.

The Arabian Sea was defined as the area with SSS > 35.3 in the northern Indian Ocean and was divided into four seasonal subzones. The area with SST > 28.3 °C (10IND-1AREX) was thought not to be affected by upwelling and was defined as a single subzone throughout the year. The area with SST < 28.3 °C was thought to be affected by lower seawater and was divided into two seasons, a SW monsoon season (Jun.–Sep.) and the remainder of the year (10IND-2ARnsw). A threshold SST of 27 °C was used to further divide the former into two regions, the northwestern Arabian Sea (10IND-2ARsw1), which is greatly affected by upwelling, and the remainder of the subzone (10IND-2ARsw2), which is only moderately affected by upwelling. We distinguished the equatorial subzone into a northern Indian Ocean area characterized by SSS < 35.3 or an equatorial area characterized by SSDH < 0.85 m (10IND-3EQ).

The rest of the southern Indian Ocean was divided into two subzones, a western subtropical subzone (10IND-4WSBT) and an eastern subtropical subzone (10IND-5ESBT) distinguished by SSDH = 0.85 m. The northern, equatorial, and Southern Indian Ocean were divided by using the SSDH minimum along the equator and the SSDH maximum along the Tropic of Capricorn.

2.5 Southern Ocean

The Southern Ocean overall is a large CO2 sink where active CO2 uptake occurs, especially in the subantarctic region during the austral summer, but the high-latitude area releases CO2 into the atmosphere during ice-free seasons (e.g., Lenton et al. 2013; Landschützer et al. 2015). The TA zones 11SOA, 12SOI, and 13SOP were first distinguished based on a nutrient-depletion line corresponding to SST = 16 °C (11SOA-1SSBT, 12SOI-1SSBT, and 13SOP-1SSBT) in a manner similar to the subtropic-like subzones. The rest of these zones and TA zone 14SSO were further divided into 0.2 m SSDH bands (11SOA-2 to 4, 12SOI-2 to 5, 13SOP-2 to 5, 14SSO-1 to 6) based on the annular variability of the carbonate system in the Southern Ocean (Takahashi et al. 2009). They were divided into austral winter (May–September) and summer (October–April) seasons based on the availability of Chl data. In the Southern Ocean, MLD was not significantly correlated with nDIC residuals (Fig. 17), and we did not use it in the MLR.

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Iida, Y., Takatani, Y., Kojima, A. et al. Global trends of ocean CO2 sink and ocean acidification: an observation-based reconstruction of surface ocean inorganic carbon variables. J Oceanogr 77, 323–358 (2021). https://doi.org/10.1007/s10872-020-00571-5

Download citation

  • Received:

  • Revised:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1007/s10872-020-00571-5

Keywords

Navigation