Development of an atmosphere-ocean coupled operational forecast model for the Maritime Continent: Part 1-Evaluation of ocean forecasts

This article describes the development and ocean forecast evaluation of an atmosphere-ocean coupled prediction 10 system for the Maritime Continent (MC) domain, which includes the eastern Indian and western Pacific Oceans. The coupled system comprises regional configurations of the atmospheric model MetUM and ocean model NEMO, at a uniform horizontal resolution of 4.5 km x 4.5 km, coupled using the OASIS3-MCT libraries. The coupled model is run as a pre-operational forecast system from 1 to 31 October 2019. Hindcast simulations performed for the period 1 January 2014 to 30 September 2019, using the stand-alone ocean configuration, provided the initial condition to the coupled ocean model. This paper details 15 the evaluations of ocean-only model hindcast and 6-day coupled ocean forecast simulations. Direct comparison of sea surface temperature (SST) and sea surface height (SSH) with analysis as well as in situ observations are performed for the ocean-only hindcast evaluation. For the evaluation of coupled ocean model, comparisons of ocean forecast for different forecast lead times with SST analysis, and in situ observations of SSH, temperature and salinity have been performed. Overall, the model forecast deviation of SST, SSH, and subsurface temperature and salinity fields relative to observation is within acceptable error limits 20 of operational forecast models. Typical runtimes of the daily forecast simulations are found to be suitable for the operational forecast applications.

Work to develop and first-hand evaluation of the WMCao, described in T18, was a preliminary step aimed to establish a highresolution atmosphere-ocean coupled model focusing on the Southeast Asia region for both operational forecast and climate 75 research applications. Since the overall objective of the development was to simulate both atmospheric and oceanic variables, the coupling has provided a better consistency between the atmospheric conditions and that of the ocean underneath rather than employing stand-alone models. The case studies conducted as part of the WMCao evaluation suggested that for an accurate prediction of weather events such as cold surges or typhoons, the zonal extent of the domain might not be sufficient.
For instance, cyclogenesis inside the South China Sea (SCS) is rare and most of the cyclones/typhoons that appear over the 80 SCS are originated in the northwestern Pacific Ocean. The north Pacific Ocean region between 100⁰E and 180⁰ is the most active tropical cyclone basin in the earth and it accounts for about one-third of the world tropical cyclone annually. Hence, rather than internal coupled dynamics, the predicted track and typhoon characteristics are dominantly driven by the Lateral Boundary Conditions (LBC) in WMCao. Similarly, the simulation of MJO and cold surge related weather parameters may also be heavily influenced by LBC. Hence, to address the issues encountered in T18 and incorporate the latest model scientific 85 developments, the present study aims to bring several key upgrades to the WMCao configuration, and test its feasibility in operational forecast application. Main updates to the coupled modelling system include extending the eastern boundary of the model domain to the west Pacific Ocean, upgrading MetUM to the latest science configuration and incorporating tide boundary forcing to the NEMO.
The study presents details of the atmosphere-ocean coupled prediction system developed for the MC and an evaluation of the 90 ocean forecast from the system using a 6-day pre-operational forecast for October 2019. Following the method employed in many earlier coupled modelling studies (e.g. Li et al., 2014;Lewis et al., 2018), the evaluation of WMCao in T18 has been performed by using short case study simulations, spanning over 5-days, of selected weather events. In the present study, instead of case studies, we assess surface and subsurface oceanic variables predicted by the coupled system across different forecast lead times. 95 The next section of this paper presents an overview of the model setup, including a brief description of the model domain, atmospheric, ocean and coupled model configurations. A brief discussion of the pre-operational forecast system setup is also presented. Section 3 provides the details of datasets used for the atmosphere/ocean model forcing and evaluation. Section 4 presents an assessment of the sea surface variables simulated by the stand-alone ocean model and both surface and subsurface ocean forecast delivered by the MC coupled model. Finally, section 5 summarizes the results obtained from the study and 100 suggests future developments.

Model domain
The model domain extends from 92⁰E to 141⁰E and 18⁰S to 24⁰N (Figure 1) on a regular latitude-longitude grid, that covers most of the tropical regions of eastern Indian and western Pacific Oceans. The deepest oceanic trench on the Earth, known as the Mariana Trench, is located in the northwestern Pacific Ocean. The crescent-shaped trench is positioned roughly between 140⁰E, 10⁰N and 150⁰E, 60⁰N (Gvirtzman and Stern, 2004). The model eastern boundary is limited to 141⁰E to avoid numerical instabilities that may arise due to steep bathymetric slopes such as the Mariana Trench. Both the atmospheric and ocean components of the coupled system are selected to have the same domain. The horizontal resolution of the MC coupled model retained to be same (4.5 km x 4.5 km) as that of the WMCao. The MC atmosphere-ocean coupled model configuration is 110 referred to as MCao in this manuscript. .

Atmospheric model
The atmospheric component of T18 has been improved to employ the SINGV v5 science configuration described in Huang et al. (2019), which is similar to the Regional Atmosphere v1 in the Tropics (RA1T) configuration of MetUM as described in Bush et al. (2020). The model is employed operationally by the Meteorological Service of Singapore since 2019 at a higher 115 resolution (1.5 km) for the region 95⁰E to 109⁰E, 6⁰S to 8⁰N and is referenced in the literature as SINGV. The key differences of MCao atmospheric model component to T18 are;  Lateral boundary conditions are provided at 3 hourly frequency from the deterministic ECMWF forecasts instead of the MetUM global deterministic model. This change has led to a significant increase in precipitation skill scores across all spatial scales and precipitation thresholds in SINGV (Huang et al., 2019); 120  The model uses a Prognostic Cloud fraction and Prognostic Condensate scheme (PC2, Wilson et al., 2008), instead of the diagnostic scheme of Smith (1990). This change helped to reduce the occurrence of spurious convection with very high rainfall rates and resulted in a better organization of convection, as shown in Dipankar et al. (manuscript submitted) The rest of the model formulations are similar to T18 and the SINGV configuration as described in Huang et al. (2019). Main 125 characteristics of the model are summarized below;  The dynamical core is the non-hydrostatic semi-Lagrangian and semi-implicit Even Newer Dynamics for the General Atmospheric Modelling of the Environment (ENDGAME, Wood et al., 2014), with an Arakawa C staggered grid.
The model time step is 120 s.  A terrain-following vertical coordinate with a resolution of 80 levels and a top lid at 38.5 km. Vertical resolution is 130 of 5 m at the boundary layer and 1.45 km below the model top, similar to the SINGV configuration.
 Boundary Layer parametrization is based on a blending between the one-dimensional scheme of Lock et al. (2001) and the three-dimensional Smagorinsky-Lilly scheme (Lilly, 1962), blending is described in Boutle et al. (2014).
 Microphysics scheme is based on Wilson and Ballard (1999) with prognostic rain formulation and improved particle size distribution for rain as in Abel and Boutle (2012).
 Radiation scheme is based on the Edwards and Slingo (1996) scheme, with six bands in the shortwave and nine bands in the longwave (Manners et al., 2011).
 The Joint UK Land Environment Simulator (JULES, Best et al., 2011) land surface scheme with 9 surface fraction types.
 The moist conservation scheme as described in Aranami et al. (2015). 140 Atmospheric component of the MCao employed in this study is referred to as MCAao hereafter.

Ocean model
A regional version of Ocean PArallelise ocean engine within the NEMO (version 3.6_stable, revision 6232, Madec et al., 2016) framework is employed as the oceanic component of the MCao. NEMO is a primitive equation, hydrostatic, Boussinesq ocean model extensively used in climate and operational forecast applications. The MCao ocean configuration shares many features 145 of its predecessor, WMCao. Hence, only key features of the NEMO and main updates of MCao configuration are discussed here.
The model horizontal grid is in orthogonal curvilinear coordinates, with Arakawa C-grid staggering. The bathymetry of MCao is based on the General bathymetric Chart of the Oceans (GEBCO2014) 30-arc second data. The model has 51 vertical levels in terrain-following coordinate system and uses the stretching function by Siddorn and Furner (2013). The stretching function 150 maintains a near-uniform surface cell thickness (< 1 m) and hence ensures the consistent exchange of air-sea fluxes over the domain, which is critical in the atmosphere-ocean coupling. Non-linear free surface following the variable volume layer formulation by Levier et al. (2007) is used for model free surface computation. The ocean model configurations used in our study have baroclinic and barotropic time steps of 120 s and 8 s, respectively.
The Generic Length Scale (GLS) turbulence model (Umlauf and Burchard, 2013) with K-ε turbulent closure scheme and the 155 stability function from Canuto et al. (2001) are used to compute the turbulent viscosities and diffusivities. Background vertical eddy viscosity and eddy diffusivity coefficients are set to a lower value of 1.2 x 10 -6 in MCao, whereas these coefficients were 1.2 x 10 -4 and 1.2 x 10 -5 , respectively, in the WMCao. Additional vertical mixing resulting from internal tide breaking is parameterized in the model as proposed by St. Laurent et al. (2002). Both energy and enstrophy conserving scheme is used for the momentum advection. For lateral tracer diffusion, the Laplacian operator along geopotential levels with a coefficient of 20 160 m 2 s -1 is used, while iso-level bilaplacian viscosity with a coefficient of -6 x 10 7 m 2 s -1 is applied for the momentum mixing.
Implicit form of non-linear parameterization with a log-layer formulation is used for the bottom drag coefficient computation.
The minimum and maximum of the drag coefficient are set to 0.0001 and 0.15, respectively.
At the lateral open ocean boundaries, the flow relaxation scheme (FRS, Davies, 1976) is applied for the tracers and baroclinic velocities, while Flather boundary condition (Flather, 1976) is used for the sea surface height (SSH) and barotropic velocities.
One of the key updates to MCao is the implementation of tide forcing at the lateral boundaries and tide potential at the ocean surface. Due to certain numerical issues, the tide related forcings are not included in the WMCao. The tidal elevations and currents from Finite Element Solutions (FES2014b) data have been used for providing the tidal harmonics at the lateral boundaries. Fifteen major tidal constituents (Q1, O1, P1, S1, K1, 2N2, Mu2, Nu2, N2, M2, L2, T2, S2, K2 and M4) are included in the boundary forcing. 170 Both coupled and uncoupled ocean model configurations are employed in the study. For uncoupled simulations, the air-sea heat fluxes are estimated using the Common Ocean-ice Reference Experiment (CORE) bulk formulae (Large and Yeager, 2004). However, a direct flux formulation is used in the coupled ocean model. Monthly runoff climatology from Dai and Trenberth (2002) and chlorophyll monthly climatology from SeaWiFS satellite observation are provided as runoff forcing and to compute light absorption coefficients, respectively, in all ocean configurations. The Red-Blue-Green (RGB) scheme is used 175 to calculate the penetration of shortwave radiation into the ocean (Lengaigne et al., 2007). Identical to WMCao, the fraction of solar radiation absorbed at the surface layer is defined to be 56% of the downward component. Mean sea level pressure (MSLP) forcing is included in the surface boundary forcing to take account of the inverse barometric effect on SSH.
The uncoupled and coupled ocean model configurations employed in this study are referred to as MCO and MCOao, respectively. 180

Coupled configuration
The exchange of fluxes between the atmosphere and ocean models is achieved through the Ocean Atmosphere Sea ice Soil coupler (version 3.3) interfaced with the Model Coupling Toolkit (OASIS3-MCT) libraries (Valcke, 2013). The Earth System Modelling Framework (ESMF) regrid tools are used to generate the interpolation weights for the remapping of exchange fields.
The coupling occurs at hourly frequency and hourly mean fields are exchanged. Since a direct flux formulation is implemented, 185 the heat fluxes computed using the Monin-Obukhov similarity theory is exchanged from the atmosphere to the ocean model.
The sea surface temperature (SST) and zonal and meridional surface current fields are sent from the ocean to the atmosphere model. The variables exchanged from atmosphere to the ocean include; non-solar heat flux, net shortwave radiation, liquid precipitation, net evaporation, and zonal and meridional wind stress. Due to numerical issues, MSLP exchange from the atmosphere is not enabled in the MCao. Instead, it is supplied from an external data source to the ocean model. 190

Model initialization and forcing
To assess the performance of the ocean model and provide initial condition to the MCOao, a 69-month hindcast simulation is performed with MCO for the period 01 January 2014 to 30 September 2019. The MCO is initialized in 1 January 2014 using temperature, salinity, zonal and meridional currents, and SSH derived from Mercator global ocean reanalysis. The lateral boundary condition for the hindcast simulation is also obtained from the same ocean reanalysis data. The daily mean of 195 temperature, salinity, baroclinic and barotropic velocities, and SSH are included in the lateral boundary forcing. Ocean surface   A schematic of the atmosphere-ocean coupled system used in the pre-operational forecast is shown in Figure 2b. In the coupled prediction system, the MCAao is initialized daily at 00Z UTC from the ECMWF IFS analysis. MCO run for the previous day 205 (T0 minus 1), forced by 6-hourly ECMWF IFS analysis at the surface and the daily mean of updated Mercator ocean forecast as the LBC, provides the initial condition to MCOao. Since it is driven by analyzed (or updated) surface (lateral) boundary conditions, the MCO provides an updated initial condition to the MCOao daily. The MCao forecast run is driven by LBC from 3-hourly ECMWF IFS forecasts in the atmosphere and daily Mercator forecasts in the ocean. Since the MSLP from MCAao is not incorporated in MCOao, 3-hourly EMCWF IFS forecast data is supplied to the model. 210

Pre-operational forecast setup
The atmosphere-ocean coupled forecast model has run as a pre-operational forecast system from 1 October 2019 to 31 October 2019 at the Cray XC-40 HPC located in the Center for Climate Research Singapore (CCRS), Singapore. The forecast system includes all necessary programs/scripts for the pre-processing of atmospheric and oceanic variables to their respective model grids. The forecast system is scheduled to initialize the forecasts daily at 1300 UTC and simulations are completed by ~1840 215 UTC. Summary of HPC resources usage and typical runtimes for daily forecast simulations are shown in Table 1. To minimize the output size, only basic oceanic and atmospheric variables are included in the output. The forecast from MCOao includes instantaneous SSH, hourly averaged sea surface temperature, sea surface salinity, and surface current velocities, and daily mean of ocean temperature, salinity and ocean currents. Further, to test the feasibility of the coupled forecast system for operational purpose, we have conducted simulations with increased computational resources. Test simulations showed that by 220 increasing the computational resources to 81 nodes (2916 cores), the total runtime has been reduced to ~140 min. This suggests a near-linear reduction in total runtime with an increase in computation nodes.

Configuration
Uncoupled ocean

Data 225
A brief description of the reanalysis, forecast and observational datasets used for the model initialization, forcing and evaluation is presented in this section.

Model initialization and forcing
ERA5 is a climate reanalysis produced by the ECMWF providing hourly estimates of many atmospheric, land and oceanic fields (Hersbach et al., 2020). Currently, it covers the period from 1979 to within 5-days of present-time and horizontal 230 resolution is approx. 30 km. The reanalysis produced using 4D-Var assimilation of ECMWF Integrated Forecast System (IFS).
ERA5 combines vast amounts of historical observations into global estimates using advanced modelling and data assimilation systems. The data is freely available through the data server https://cds.climate.copernicus.eu/.
ECMWF IFS is a global weather prediction system comprising a spectral atmospheric model, ocean wave model, ocean model and land surface model coupled to a 4D-Var data assimilation system. IFS medium-range weather forecasts are available up 235 to 10 days at a horizontal resolution of 0.1⁰. In addition, the atmospheric analysis fields are provided four times daily for the forecast base time 00, 06, 12 and 18 UTC. The data is available to registered users from https://www.ecmwf.int/en/forecasts/datasets/ Mercator global ocean reanalysis/forecast provides oceanic variables with 1/12⁰ horizontal resolution. The system uses NEMO v3.1 with 50 vertical z-levels ranging from zero to 5500 m and forced by the ECMWF IFS meteorological variables. The 240 assimilation/forecast product includes the daily mean of temperature, salinity, currents from top to bottom over the global ocean and SSH. The data is freely available from https://marine.copernicus.eu/.
The tidal heights and currents computed from the global tide model Finite Element Solution (FES2014b) is used as the tidal forcing in the model. FES2014 is based on the resolution of the shallow water hydrodynamic equations (T-UGO model) in a spectral configuration and using a global finite element mesh with increasing resolution in coastal and shallow waters regions 245 (Lyard et al., 2006). The database is distributed on a global 1/16° x 1/16° grid. Data is produced by assimilating long-term altimetry data (Topex/ Poseidon, Jason-1, Jason-2, TPN-J1N, and ERS-1, ERS-2, ENVISAT) and tidal gauges through an improved representer assimilation method. Tidal heights and currents of 32 tidal constituents are available. The data is freely available through http://www.aviso.altimetry.fr/en/data/products/auxiliary-products/global-tide-fes.html.

Model evaluation 250
The CORIOLIS data service provides quality-controlled is situ data in real-time and delayed mode over the global ocean. The data include temperature and salinity profiles and time series from profiling floats, Expendable Bathy Thermograph's (XBT), thermo-salinographs (TSG), and drifting buoys. The data is freely available from http://www.coriolis.eu.org/Data-Products/.
The Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) SST is produced daily on an operational basis at the UK Met Office using optimal interpolation on a global 0.054⁰ x 0.054⁰ grid. The product assimilates satellite data including 255 advanced very-high resolution radiometer, Spinning Enhanced Visible and Infrared imager, Geostationary Operational Environmental Satellite Imager, Infrared Atmospheric Sounding Interferometer, Tropical Rainfall Measuring Mission Microwave imager and in situ data from ships, drifting and moored buoys (Donlon et al., 2012). SST data at every grid point is accompanied by an uncertainty estimate, known as an analysis error, and an optimal interpolation approach is employed to produce this estimate. The data is freely available from https://marine.copernicus.eu/. 260 The University of Hawaii Sea Level Center (UHSLC) offers quality controlled tide gauge (TG) sea level observations over the global ocean as fast delivery (FD, 1-2 months delay) and research quality (RQ, 1-2 year delay) data at hourly and daily resolution (Caldwell et al., 2015). The data is freely available from http://uhslc.soest.hawaii.edu/data/.

Ocean Hindcast
An overall assessment of the ocean model hindcast simulation is carried out to understand the realism of the ocean initial  Table 2. SST Bias, RMSD and correlation coefficient statistics computed using modelled SST and OSTIA are shown in the Table. The SST bias is within +0.2 ⁰C for about 76% and within +0.5 ⁰C for about 98% of the analysis domain.
Largest SST cold bias is seen in the Andaman Sea region. Meanwhile, most of the South China Sea (SCS), equatorial west 295 Pacific Ocean and Australian coasts of the Timor Sea show a positive SST bias. Negative SST bias of about -0.25 ⁰C is observed in the ASMS region, while positive bias over 0.25 ⁰C is confined to the SSCS and GoT sub-regions. Rather than appearing as a basin-wide feature, higher positive biases appear as small circular patches in the northern SCS region that represents the likely existence of cyclonic eddies over this region.   The RMSD between model and OSTIA is less than 0.5 ⁰C for about 97% of the analysis-domain (Figure 4b). Small patches of higher RMSD (> 0.7 ⁰C) are mostly seen along the coastal regions. The RMSD minimum (0.25 ⁰C) and maximum (0.53 ⁰C) are observed over the BS and ASMS sub-regions, respectively (Table 2). Correlation between the model SST hindcast and 305 OSTIA is above 99.9% confidence level over the analysis-domain (Figure 4c). Over 88% of the domain displays a correlation higher than 0.8. Relatively low correlation is seen over the middle of Malacca Strait, Makassar Strait and equatorial Pacific Ocean regions. In sub-region spatial average, the lowest (0.8) and highest (0.96) correlations are seen over the SuCeS and RSCS regions, respectively (Table 2). Time series of spatially averaged SST difference between model and OSTIA is shown in Figure 4d. Consistent with our earlier analyses, relatively low SST difference depicts a good agreement between the MCO 310 SST hindcast and OSTIA analysis.  Temperature observation at 1 m depth is taken as SST from the moored buoys while temperature averaged over upper 1 m is indicated as the model SST at these locations. In general, a good agreement is found between the model and observations at 320 both mooring locations. Both the seasonal and intra-seasonal SST variability is reasonably well reproduced by the model. SST bias, RMSD and correlation between the model and observation are 0.17 ⁰C, 0.29 ⁰C and 0.94, respectively, for M1 and 0.12 ⁰C, 0.41 ⁰C and 0.92, respectively, for M2 locations. The standard deviation (SD) of SST at M1 and M2 locations are 0.94 ⁰C and 0.99 ⁰C, respectively, and the RMSD is relatively smaller than the SD at both locations.

Sea surface height 325
Daily mean SSH observation from 20 tide-gauge stations distributed across the domain and MCO simulated SSH interpolated to the location of these observations have been used for the hindcast evaluation. SSH bias, RMSD and correlation coefficient statistics between model and SSH observations are given in Table 3   Time series of daily mean SSH from model and observations for randomly selected stations are plotted in Figure 6. Where the stations Sibolga and Prigi are located in the eastern tropical Indian Ocean, and Currimao Ilocos Norte and Vung Tau are located in the SCS. Generally, the model simulated SSH follows the observation and shows good agreement with it. A few sharp peaks in the tide-gauge observation are found to be absent in the model simulation (e.g. Figure 6c). Most of the tide gauge stations 340 are located adjacent to the coast and our current model resolution is not enough to resolve the coastline at very fine scales.
Since the model SSH is interpolated to the tide gauge location, local scale SSH variations may not be captured in the model simulation. This may be one possible reason for the discrepancy between model and observation.

Ocean Forecasts
Results from the analysis of MCOao forecast simulations for October 2020 are presented here. Since the system delivers a 6-350 day forecast, the analysis period extends from 1 October to 5 November 2020. Daily files are produced at different forecast lead times, T+0 to T+24 (fsct_day1), T+24 to T+48 (fcst_day2), T+48 to T+72 (fcst_day3), T+72 to T+96 (fcst_day4), T+96 to T+120 (fcst_day5) and T+120 to T+144 (fcst_day6) for the following analyses. Comparisons of coupled ocean forecast for different forecast lead times with OSTIA SST and in situ observations such as temperature from RAMA moored buoys, TSG and XBT profiles, temperature and salinity from Conductivity Temperature Depth (CTD) profiles, and SSH from tide-gauges 355 have been performed.

Sea surface temperature
Time series of daily mean SST averaged over the sub-regions from OSTIA analysis and different forecast lead times are plotted in Figure 7. Statistics of SST bias, RMSD and correlation coefficient between model and OSTIA for the October forecast run are listed in Table 4. The forecasted SST over most of the sub-regions is within the error standard deviation of the OSTIA 360 analysis, which is indicated by shading in Figure 7. Excluding the ASMS, all other sub-regions exhibit a warm SST bias with the largest values over the SSCS and GoT. The RMSD is less than 0.5 ⁰C over most of the sub-regions during the analysis period. Over ASMS sub-region, both the cold bias and RMSD increases with the forecast lead time and shows the highest RMSD (0.49 ⁰C) on fsct_day6. The largest SST RMSD (0.56 ⁰C) and bias (0.49 ⁰C) over the entire sub-regions is observed over the GoT on 1-day forecast lead time. Generally, the forecasted SST tends to be cooler with an increase in forecast lead 365 time, denoting a lower warm bias and RMSD relative to fcst_day1. Interestingly, inconsistent with the improvements in SST bias and RMSD, the correlation between forecasted and OSTIA time series considerably decreases with higher forecast lead times (Table 4).
SST correlation is above 95% confidence level (r > 0.365) over entire analysis-domain during fcst_day1. In general, SST correlation is higher than 99% confidence level (r > 0.46) over 60% of the sub-regions during higher forecast lead times as 370 well. The sub-regions including ASMS, GoT, TWPO and SuCeS are noted by lower correlation significance level with an increase in forecast lead time. Overall, the SST correlation over the analysis-domain is above 99% confidence level across all forecast lead times while bias and RMSD are less than 0.19 ⁰C and 0.35 ⁰C, respectively.   Despite this, the correlation between model forecasts and observation depicts a considerable decrease at higher forecast lead times. A cold SST bias of about -0.14 ⁰C to -0.17 ⁰C is noted at location M2. The RMSD at M2 is relatively less, with a maximum of 0.18 ⁰C across all forecast lead times while compared to M1. The SST correlation at M2 is above 95% confidence level (r > 0.61, df > 9) during all forecast lead times.    Temperature and salinity data available from the CORIOLIS data portal is employed for the MCOao forecast evaluation of subsurface fields. Argo and XBT profiles, moored buoys, and TSG observations are used for the analysis. We only consider mooring and TSG observations with a minimum of 10 days for the analysis. TSG observation selected for the analysis is located in the Timor-Arafura Sea and the track has a direction of motion towards the west (Figure 3). Continuous observation 400 available at 6.4 m depth is compared with the model forecasts. Currently, the forecast system produces only daily averaged subsurface variables. Hence, the instantaneous TSG temperature observations at 1200 UTC are compared with the daily averaged model temperature. Time series of TSG temperature observation and model forecast for different forecast lead times is shown in Figure 8c. Though the model temperature shows a negative bias, daily temperature variation is reasonably well predicted by the model. Which is supported by high correlation, above 99.9% confidence level, between the model forecast 405 and observation (Table 5, item no. 3). Largest temperature bias and RMSD are -0.47 ⁰C and 0.68 ⁰C, respectively, and significant variations in bias, RMSD and correlation statistics are not noticeable.
The moored buoy observations provide a unique opportunity to assess the model simulations both in the surface and subsurface of the ocean. Since the variables of interest in the present study are temperature and salinity, we have tried to compare the model forecast of temperature and salinity with buoy observation available at the locations M1 and M2. However, due to data 410 gaps and shorter time series, salinity observation from both moorings and temperature from M2 mooring are not included in the analysis. Depth-time plots of temperature from the model for forecast lead times 1 day (fcst_day1), 3 days (fcst_day3) and 5 days (fcst_day5), and moored observation, M1, are shown in Figure 9. Statistics of temperature bias, RMSD and correlation coefficient between model forecast and observation are given in Table 5. The depth of the upper ocean isothermal/mixed layer and its shoaling in late October are well simulated by the model. A significant difference between the model forecast and 415 observation is seen in the region below the mixed layer. The model simulation is unable to reproduce the sharp temperature stratification/cooling in the thermocline regions while this feature is well evident in the observation. This leads to relatively larger discrepancy between the model and observation in the subsurface region roughly between 100 m to 250 m depths. The usage of daily averaged temperature rather than instantaneous profile or higher vertical mixing in the model may be one of the possible reasons for this discrepancy. Larger temperature difference at the thermocline region has led to a warm temperature 420 bias in the model forecast (Table 5). Maximum temperature bias and RMSD are 1.95 ⁰C and 2.96 ⁰C, respectively. Meanwhile, the correlation between the model forecast and observation is above 99% confidence level (r > 0.47) across all forecast lead times.
Argo and XBT profiles available for the period 1 October to 5 November 2019 are compared with model forecast to derive the RMSD statistics for temperature and salinity. Since no temperature or salinity profiles are available in the SCS (figure not  425 shown, data distribution can be viewed from http://www.coriolis.eu.org/Data-Products/Data-Delivery/Data-selection), the analysis mainly demonstrates the model performance in the domain excluding the SCS region. As observed in the M1 mooring location, warm biases with varying magnitude are seen, in the thermocline region, across the analysis-domain (figures not shown). Considering the depth range where this bias exists, the vertical mixing parameterization may have a strong influence in modifying the thermal stratification than the penetrative short wave forcing. Statistics of RMSD for ocean temperature and 430 salinity relative to all profile observations is given in Table 6. RMSD of individual profiles are first computed and then rms value of the computed RMSD is derived. RMSD across all forecast lead times along with the number of profiles analyzed are listed. For both temperature and salinity, the RMSD remains fairly similar during the entire analysis period. Over the analysis-domain and across all forecast lead times, the maximum RMSD for temperature and salinity are 1.41 ⁰C and 0.14 psu, respectively. Overall, the model forecast deviation relative to observation is within acceptable error limits of operational 435 forecast models (e.g. Zhang et al., 2010;Yang et al., 2016).

Sea surface height
The same set of tide-gauge stations used for the MCO hindcast validation has been employed for MCOao SSH forecast evaluation. Due to irregularities in time series, the Currimao Ilocos Norte station is not included in the analysis. Hourly 445 instantaneous SSH data from the model forecast and tide-gauge observation is used for the statistical analysis. Summary of SSH RMSD and bias statistics relative to the observations are listed in Table 7. Time series of hourly instantaneous SSH at randomly selected tide-gauge stations (Sibolga, Prigi and Vung Tau) and MCOao forecast (fcst_day1) are plotted in Figure 10.
Model SSH bias is within +0.10 m at all tide-gauge stations and forecast lead times. The SSH bias is within + 0.05 m for 14 of total 19 tide-gauge stations across all forecast lead times. Since the tide-gauges are mostly located near to the coast, SSH 450 variability shorter than intra-seasonal timescale may be largely driven by the tidal forcing. Eventually, in the SSH bias, there will be an offset between high and low tidal peaks. Hence, RMSD will give a better representation of model accuracy in tidal  Mean RMSD and Bias 0.14 0.14 0.14 0.14 0.14 0.14 0.00 0.00 0.00 0.00 0.00 0.00

Summary and future developments
The Maritime Continent has a profound influence on the global climate system because of its complex topography and unique 465 geographic location within the tropical Indo-Pacific warm pool. The MC region is characterized by strong atmosphere-ocean coupled processes across multiple timescales. A convective-scale/eddy-resolving atmosphere coupled modelling system for the western MC, described in T18, was able to improve the simulation of a cold surge event, the intensity of Typhoon Sarika and its atmosphere-ocean interactions. Several upgrades has been added to the T18 model for a future operational implementation, such as extending the eastern boundary of the model domain to the west Pacific Ocean, improved science 470 configuration of the atmospheric model, MetUM, and incorporation of tidal boundary forcing to the ocean model, NEMO.
Furthermore, the coupled model's feasibility to use as an operational forecast system is also being tested. Typical runtimes of the daily forecast simulations in our experiments are found to be suitable for the operational forecast applications.
The MCao coupled prediction system has run as a pre-operational forecast system from 1 to 31 October 2019. Hindcast simulations performed for the period 1 January 2014 to 30 September 2019, using the uncoupled ocean model MCO, has 475 provided the initial condition to the MCOao. The paper present details of atmosphere-ocean coupled prediction system developed for the MC, and evaluations of ocean-only model hindcast and 6-day ocean forecast simulation performed using the coupled system. the OSTIA. Though the model forecast exhibits a warm SST bias, the RMSD is less than 0.45 ⁰C over most of the sub-regions during the analysis period. Generally, the forecasted SST tends to be cooler with an increase in forecast lead time, denoting a lower warm bias and RMSD relative to fcst_day1. Overall, the SST correlation over the analysis-domain is above 99% confidence level across all forecast lead times while the bias and RMSD are less than 0.19 ⁰C and 0.35 ⁰C, respectively.
The diurnal variability of SST at the RAMA moored buoy locations M1 and M2 are reasonably well reproduced by the model 500 across all forecast lead times. SST bias and RMSD at M1 are less than 0.07 ⁰C and 0.20 ⁰C, respectively, and remains fairly constant across all forecast lead times. Meanwhile, the RMSD at M2 is relatively less, with a maximum of 0.18 ⁰C across all forecast lead times. The depth of the upper ocean isothermal/mixed layer at the M1 is well forecasted by the model. To understand the model skill in predicting subsurface temperature and salinity, in situ profile observations are compared with the model forecast for different forecast lead times. For both temperature and salinity, the RMSD remains fairly constant during 505 the entire analysis period.
Comparison of model forecasted SSH shows good agreement with the tide-gauge observations. About 75% of the stations show bias within + 0.05 m at all forecast lead times. No significant variation in model forecast accuracy or RMSD is seen with the increase in forecast lead time. About 73% of total stations show RMSD less than 0.10 m for all forecast lead times. Overall, the model forecast deviation of SST, SSH, and subsurface temperature and salinity fields relative to observation is within 510 acceptable error limits of operational forecast models (e.g. Zhang et al., 2010;Yang et al., 2016).
Analysis of subsurface fields revealed that significant model temperature biases exist in the region below the mixed layer. The representation of stratification in the thermocline region is relatively weak in the model forecast than the observations. This subsurface temperature bias remains across the model domain with varying magnitudes. Similar temperature biases are earlier reported in simulations with identical model configurations (e.g. Graham et al., 2018). Further modelling works are needed to 515 improve the thermal stratification through fine-tuning the mixing coefficients or modifying the vertical mixing parameterization. Further analysis of model forecast fields using longer forecast simulations with increased observations will be undertaken to assess the model predictability across different seasons and during typical weather events such as cold surge, typhoon, MJO.
In addition, the impact of coupling on the forecast will be investigated by performing simultaneous stand-alone ocean model 520 forecast simulations. Further, works are ongoing to understand the forecast skill of MCAao and results from the analysis will be presented as research publications.
The evaluation of sea surface salinity and ocean current fields are not included in the present study mainly due to the lack of in situ and satellite observations. Especially, the South China Sea remains as a data-sparse region in terms of the subsurface observations that highlights the necessity of coordinated efforts from the scientific community to fill these spatial data gaps. 525 In our analyses, the RMSD and positive or negative biases of the ocean forecasts are generally comparable to those observed from the hindcast statistics. This suggests that up to a certain extent, the model forecast deviation is inherited from the MCA hindcast or the MCAao initial condition. The dependency of model forecast quality on the initial state is well established in the numerical weather prediction studies. In particular, over the tropical oceans, the initialization of the ocean state is an important element of the forecast systems. The data assimilation techniques help to acquire an improved estimate of the ocean state by 530 combining the model simulated fields and observations (King et al., 2018). Both the uncoupled and coupled ocean configurations used in our study are free-running models with no restoring or relaxation to the real world. Hence, to provide a better initial condition to MCOao, implementation of data assimilation capability to MCO will have a key priority in our future developments. Also, earlier studies have shown that including wave-induced mixing in ocean circulation model yields a better representation of the upper ocean temperature (Lewis et al., 2019b). Thus, work towards the development of a two-way 535 atmosphere-ocean-wave coupled system will be undertaken in the future. https://doi.org/10.5194/gmd-2020-326 Preprint. Discussion started: 19 October 2020 c Author(s) 2020. CC BY 4.0 License.

Model code and data availability:
Due to intellectual property right restrictions, the coupled model system code cannot be provided directly. However, full source codes, scripts and related forcing datasets used in the study can be made available to the Editor for review. The Met Office Unified Model (MetUM) is available for use under licence from UK Met Office via a shared MetUM code repository, which 540 can be accessed via https://code.metoffice.gov.uk/trac/um/wiki. The NEMO vn3.6 model code is freely available from the NEMO website (https://www.nemo-ocean.eu). The OASIS3-MCT coupler is disseminated to registered users as free software from https://verc.enes.org/oasis. All model outputs analyzed in the manuscript can be made available upon contacting the authors.