## 2 Methodology  
The multi-scale ocean model used in this study is an improved version of the original Maritime Environment Committee (MEC) ocean model developed by the Japan Society of Naval Architecture and Ocean Engineers (JASNAOE). In the MEC model, hydrostatic approximation and non-hydrostatic(i.e., full-3D) models are applied for the meso- and small-scale domains, respectively. In the mesoscale domain, tidal flow is generated under hydrostatic (pressure) approximation to reduce the computational time. The applicability and accuracy of the hydrostatic model of the MEC model can be found in the research of [PERSON] et al. [11], [PERSON] et al. [14], and so on. A full-3D small-scale model is necessary when the vertical flow component cannot be ignored. The spatial connection between the two models is such that the full-3D model domain matches one grid-column of the hydrostatic model. [PERSON] et al. [11] modified the MEC model by adopting a two-phase model for the small-scale domain, where the continuous liquid phase (seawater) and dispersed gas/liquid phase (individual CO\({}_{2}\) bubble/droplet) are solved using Eulerian and Lagrangian methods, respectively. This model, hereafter termed the \"MEC-CO\({}_{2}\) model\", was used in the present study. The governing equations of the mesoscale model adopting the hydrostatic approximation are as follows. \[\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{ \partial z}=0 \tag{1}\]  
\[\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u }{\partial y}+w\frac{\partial u}{\partial z}=-\frac{1}{\rho}\frac{\partial p }{\partial x}+fv+A_{M}\bigg{(}\frac{\partial^{2}u}{\partial x^{2}}+\frac{ \partial^{2}u}{\partial y^{2}}\bigg{)}+\frac{\partial}{\partial z}\bigg{(}K_ {M}\frac{\partial u}{\partial z}\bigg{)} \tag{2}\]  
\[\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v }{\partial y}+w\frac{\partial v}{\partial z}=-\frac{1}{\rho}\frac{\partial p }{\partial y}-fu+A_{M}\bigg{(}\frac{\partial^{2}v}{\partial x^{2}}+\frac{ \partial^{2}v}{\partial y^{2}}\bigg{)}+\frac{\partial}{\partial z}\bigg{(}K_ {M}\frac{\partial v}{\partial z}\bigg{)} \tag{3}\]  
\[0=-\frac{\partial p}{\partial z}-\rho g \tag{4}\]  
where \(u,\ v,\) and \(w\) are the velocity components in the \(x,\ y,\) and \(z\) directions, respectively. \(p,\)\(\rho,\) g, and f are the pressure, density of seawater, gravity acceleration, and Coriolis force coefficient, respectively. \(A_{M}\) and \(K_{M}\) are the horizontal and vertical eddy viscosities, respectively. The transport equations of scalar properties \(\phi\) (i.e., temperature (\(T\)), salinity (\(S\)), and DCO\({}_{2}\) (\(C\))) are as follows:  
\[\frac{\partial\phi}{\partial t}=A_{D}\bigg{(}\frac{\partial^{2}\phi}{ \partial x^{2}}+\frac{\partial^{2}\phi}{\partial y^{2}}\bigg{)}+\frac{ \partial}{\partial z}\bigg{(}K_{D}\frac{\partial\phi}{\partial z}\bigg{)} \tag{5}\]  
where \(A_{D}\) and \(K_{D}\) are the horizontal and vertical eddy diffusivities, respectively. As mentioned earlier, the Eulerian-Lagrangian two-phase model is applied for small-scale full-3D regions. The continuity and Navier-Stokes equations for the continuous phase are as follows:  
\[\frac{\partial}{\partial t}(\alpha_{c}\rho_{c})+\
abla\cdot(\alpha_{c}\rho_{c }\mathbf{u}_{c})=\frac{1}{V_{cell}}\sum_{n=1}^{n_{d}}\Gamma \tag{6}\]  
\[\frac{\partial}{\partial t}(\alpha_{c}\rho_{c}\mathbf{u}_{c})+\
abla\cdot(\alpha_ {c}\rho_{c}\mathbf{u}_{c}\otimes\mathbf{u}_{c})+\frac{1}{V_{cell}}\sum_{n=1}^{n_{d}} \frac{\text{D}_{\text{d}}}{\text{Dt}}(\rho_{d}V_{d}\mathbf{u}_{d})=-\
abla P+ \
abla\cdot[\alpha_{c}\rho_{c}(v_{c}+v_{t})\mathbf{d}_{c}] \tag{7}\]  
where,  
\[\alpha_{d}+\alpha_{c}=1 \tag{8}\]  
\[\alpha_{d}=\frac{1}{V_{cell}}\sum_{n=1}^{n_{d}}V_{d} \tag{9}\]  
\[\mathbf{d}_{c}\mathbf{=}\
abla\mathbf{u}_{c}+(\
abla\mathbf{u}_{c})^{T} \tag{10}\]  
where \(\alpha\) is the volume fraction and \(\mathbf{u}\) is the velocity vector.