This paper takes the POM with generalized coordinate system as a representative example, without loss of generality. The vertical distribution calculation scheme in this paper is not dependent on a specific type of vertical coordinate. The Cartesian coordinate system ( x , y , z , t )   can be transformed into the vertical generalized coordinate (s-coordinate) system ( x * , y * , k , t * ) through Equation 1.

where s ( x * , y * , k , t * ) = σ ( k * ) ( H ( x * , y * ) + η ( x * , y * , t ) ) , k is the vertical k th layer, 1 ≤ k ≤ k b , k b represents the lowest layer, and σ ( 1 ) = 0 , σ ( k b ) = 1 . H ( x * , y * ) can either represent a constant or serve as a function. η ( x , y , t ) is the surface elevation. H ( x , y ) is the bottom topography.

Considering the salinity function S ( x , y , z , t ) = S * ( x * , y * , k , t * ) , Equation 2 can be derived.

Taking the partial derivative of Equation 1 with respect to SS, we obtain Equation 3.

Substituting Equation 3 into Equation 2, we obtain Equation 4.

Incorporating the aforementioned transformation process into the ocean primitive equations, ∂ S ∂ t + V → . ▿ S = k s Δ S , we can derive the basic equations used in traditional numerical models. Within those equations, the equation for salinity is represented by Equation 5:

where ω = W − ( s x + η x ) U − ( s y + η y ) V − ( s t + η t ) , K H is vertical diffusivity, q S x = A H ∂ S / ∂ x and q S y = A H ∂ S / ∂ y represents the horizontal diffusion term for salinity, A H represents the horizontal eddy viscosity mixing coefficient.

In the POMgcs, a leapfrog scheme is adopted for the solution. The computation process is divided into two steps. In the first step, the advection term and horizontal diffusion term are solved explicitly, as shown in Equation 6. In the second step, the vertical diffusion term is solved implicitly, as demonstrated in Equation 7:

Where n represents n th time step, S ˜ refers to the intermediate variable in the leapfrog scheme., a d v ( s k S ) = ∂ U s k S ∂ x + ∂ V s k S ∂ y + ∂ ω S ∂ k and d i f ( s k S ) = ( ∂ s k q S x ∂ x + ∂ s k q S y ∂ y ) represents the advection and the diffusion term of salinity respectively. The numerical discrete scheme of Equation 6 in the vertical direction can obtain Equation 8:

where k is the vertical k th layer, σ represents the horizontal area of the current layer, and z represents the thickness of the current layer. Equation 8 provides a detailed description of Equation 7. The left-hand side of the equation represents the accumulation of the product of salinity and volume across various vertical layers during an intermediate process, which equates to the total salt content of the seawater in that process. The right-hand side expresses the cumulative total salt content in each layer over time steps and the cumulative total salt content introduced (or removed) due to advection and horizontal diffusion at the time step. However, in the actual programming implementation process, to further simplify the equation, a default assumption is made. Instead of the cumulative total salt content being equal for each layer, it is assumed that the corresponding salt content in each layer is equal. Under this assumption, it is implicitly accepted that the positions of the vertical layers do not change over time. This leads to Equation 9:

Up to this point, the entire process described above pertains to the traditional method for calculating vertical salinity. It is at this juncture that default assumptions lead to artificial error being introduced. If there is no introduction of surface freshwater, this equation will not yield artificial errors. However, in the past, there have been two primary issues when dealing with this equation. One is that after introducing virtual salt flux at the boundaries to satisfy the solving conditions of the equation, the conservation of salt will be destroyed. The second is that the new stratification has deviated from the original physical positions of the previous time step. Consequently, the model does not accurately represent salinity conservation and runoff can instantaneously impact the deeper ocean. This is not natural.

In the ocean model with a free surface, the vertical layering algorithm, z = η ( x * , y * , t ) + s ( x * , y * , k , t * ) , includes the undulation of the sea surface. There are two causes for the fluctuation of the sea surface: one is due to the internal dynamic mechanisms of the ocean, that is, the convergence and divergence of seawater; the other is caused by the increase or decrease in freshwater, such as precipitation or evaporation. The assumption in Equation 9 does not introduce errors due to the former, because Equations 6 and 9 already encompass the flow across adjacent volume elements, thus ensuring the conservation of salinity in corresponding vertical layers at different time steps. However, Equation 9 fails to reflect the latter impact, namely the effects of freshwater changes. At the surface layer, we usually handle the influx and outflux of freshwater by modifying the salinity or using a virtual salt flux approach to maintain salinity conservation. But the impact of freshwater involves not only changes in physical quantities such as salinity and heat but also a shift in the position of the vertical layers. This leads to a mismatch between the actual flow positions of the seawater and the positions in the model. The former change is based on the expansion and contraction of volume elements under the premise of mass conservation, whereas the latter is the physical displacement of volume elements, which is an artificially introduced error. Consequently, this results in non-conservation of salinity and allows instantaneous and direct deep-sea impacts from phenomena like precipitation, evaporation, and runoff, which is not consistent with real-world conditions.

We split η into two parts, η * and η f , representing the sea surface height changes caused by the ocean dynamic processes and the changes due to the influx and outflux of freshwater, respectively, that is, η = η * + η f . Thus, we can derive Equation 10:

Let s d = σ ( H + η d ) , s f = σ η f . Reinserting this into Equation 1 and deriving anew, Equation 5 can be transformed into Equation 11:

where

Combining Equations 11 and 12, we can derive Equation 13. After rearranging terms, Equation 14 is obtained.

Substituting Equation 14 into the leapfrog scheme and expanding, we obtain Equation 15:

If we proceed according to Equation 9, then we can obtain Equation 16:

The influence of freshwater on marine dynamic processes is expected to be confined to the sea surface. However, Equation 16 suggests that the impact of freshwater permeates into every vertical layer, as illustrated in Figure 1 .

In the Figure 1 , from T 0 to T 1 , precipitation exceeds evaporation, leading to an increase in freshwater, and thus the actual dynamic process-induced sea surface fluctuations are lower than the traditional model stratification. Apart from the top layer, the yellow areas at the bottom should not represent freshwater with zero salinity but should be seawater. From T 1 to T 2 , precipitation is weaker than evaporation, resulting in a decrease in freshwater, so the sea surface fluctuations caused by the real dynamic processes are higher than the traditional stratification. The salt remaining after evaporation in the topmost yellow area should accumulate in the sub-surface layer. The salinity in the yellow areas of the remaining layers should not change. Therefore, we need to readjust the salinity distribution across different layers based on s * and s f to recalibrate the vertical salinity profile. Therefore, we propose a natural vertical distribution calculation scheme, as Equation 17, and then its vertical expansion yields Equations 18, 19.

where W represents the weight coefficient. The proposed scheme does not alter the conservation relationship of salinity but changes the vertical distribution of salinity. The modified salinity is then reapplied to Equation 9.

To adjust the salinity distribution after re-stratification, the weight parameter W needs to be calculated based on the real position. An indicator function can be given as Equation 20. It has been applied in Equation 21.

In order to simplify the formula for ease of writing and understanding, we have made the following simplifications, Equation 21. This simplification of Equation 21, makes the equation more concise and aesthetically pleasing, as well as easier to program.

where Z k = η + s ( k ) , Z k * = η * + s * ( k ) . Ultimately, we can derive a unified formula for the weight function.

Where d z k = s ( k ) − s ( k + 1 ) , d z k * = s * ( k ) − s * ( k + 1 ) .

Because the decomposition in η also leads to changes in other equations, it distinguishes itself from the virtual salt flux and natural boundary condition.

We treated the variation of freshwater in the surface layer as part of the ocean, encapsulated within η . Referring to previous studies (Beron-Vera et al., 1999; Warren, 2009; Tseng et al., 2016; Verri et al., 2020), for the Boussinesq continuity and conservation equations of the salt, Equations 23 and 24, their boundary conditions at the surface layer z = η are respectively Equations 25 and 26.

where, F → c H and F S V respectively represents the horizontal and vertical diffusion term of salt. This implies that salt cannot undergo any form of advection or diffusion with the atmosphere through the air-sea interface.

For the virtual salt flux, its boundary conditions we specified in our experiment are satisfied by Equations 27 and 28.

where q w is a net addition of water to the ocean ( P − E ).

For the natural boundary condition, its boundary conditions are satisfied by Equations 29 and 30.

The natural boundary condition preserves the vertical velocity in the surface, maintaining the driving effect of freshwater on seawater. The vertical diffusion term ensures the salinity balance of the water body. Our method has a similar form at z = η * . Substituting η = η * + η f into Equation 27:

We do not consider the horizontal movement of freshwater in the surface layer, u → · ∇ η f = 0 . Because of ∂ η f ∂ t = − q w , the Equations 31 and 28 can be transformed to Equations 32 and 33.

Thus, the scheme we proposed can still preserve the driving process of freshwater on the ocean.

Our method is simpler and more user-friendly compared to approaches like VLR, with less changes to the original ocean model. These modifications are straightforward to implement using the equations provided in our manuscript. For instance, considering the POM, Figure 2 depicts the flow diagram of the POM. The modifications we introduced are highlighted in red boxes within this diagram. One modification pertains to the computation of surface elevation adjustments due to freshwater fluxes in the external mode. The other critical modification involves incorporating an NVCDS computational module after the calculations for convection and horizontal diffusion, but prior to vertical diffusion. This shows that our method does not require altering the original structure of the model.

Figure 3 showcases the procedure flowchart for the NVDCS module. This module necessitates four input variables: the vertical distribution from both the preceding and current time steps, the increase in surface elevation due to freshwater fluxes, and the salinity values resulting from convection and horizontal diffusion computations. The first three variables are used to calculate weights according to Equation 22 in the manuscript, and then combined with the fourth input to compute the adjusted salinity using Equation 19. This adjusted salinity is subsequently employed in the ensuing step for vertical diffusion calculations. As illustrated, the computational sequence does not necessitate any additional reference frame conversions or interpolations, ensuring its modular integration. Furthermore, the parameters required as inputs are inherent to most existing ocean models, ensuring direct applicability.

The results presented in Figure 5 show the temporal evolution of the average salinity under different experimental conditions. It is evident that, in Ctrl experiment without implementing no salinity conservation mechanism, the salinity increases continuously over time, in sharp contrast to the other two experiments. This salinity imbalance has the potential to affect the model simulation in the long term. In contrast, the NoFlux and NVDCS experiments, which respectively do not use any salinity fluxes and use the NVDCS format to conserve salinity, show a remarkable stability of the average salinity over time. The maximum difference in the average salinity for the Ctrl experiment is 0.0017062 p s u , while the maximum difference for the other two experiments is only 4 × 10 − 13   p s u , indicating that these conservation mechanisms have a significant impact on the model simulation.

The experimental area includes three points: P1 (5°N, 30°E), P2 (31°N, 30°E), and P3 (55°N, 30°E), as well as two cross sections, ab and cd, which span the longitude and latitude of the simulated area. The location of these points is illustrated in the right panel of Figure 4 and their locations on the horizontal precipitation distribution map are shown in the same figure. These three points represent regions of high evaporation, regions where evaporation and precipitation are balanced, and regions of high precipitation, respectively. Figures 6 – 8 show the vertical salinity changes at these three points over time. It is clear that both experimental Ctrl and NVDCS have similar processes of salinity propagation, but experimental NoFlux completely loses this process. In the panels b of those figures, the early stages of precipitation are sufficiently mixed with the bottom seawater, so that it has no stratification at all in the vertical direction. Obviously, the experimental NoFlux does not use salt flux, and the salinity caused by precipitation directly affected the deep sea. In the experiment, NVDCS simulates the vertical distribution of salinity well, and obtains a salinity distribution similar to that of the virtual salt flux Ctrl.

In the NoFlux experiment, salinity changes are mainly caused by the volume change of seawater as there is no salt flux. This volume change of seawater has a direct effect on the deep sea due to coordinate stratification. As a result, there is no mixing process in the middle, and no vertical distribution occurs ( Figures 6B, E ). In the Ctrl experiment, the salinity change is caused by two factors, the change in seawater volume, and the input and output effects of salt. The former has a smaller effect, while the latter is the main contributor. Therefore, there will be significant vertical distribution ( Figures 6A, D ). However, the change in seawater volume can also result in lower salinity of deep seawater in freshwater regions during the early stages of model operation, and higher salinity of deep seawater in high evaporation regions. As the model continues to operate, and salinity continues to increase, ultimately, both precipitation and evaporation regions will generally have higher salinity in the deep sea.

In the vicinity of Point P1, evaporation is the dominant process with a higher rate than precipitation, resulting in a decrease of freshwater at the sea surface. Figure 6 depicts the distribution of salinity under different numerical experiments. Experiment Ctrl and NVDCS demonstrate a clear formation of high salinity regions at the sea surface, without significant impact on the salinity of the deep seawater. However, Experiment NoFlux does not exhibit any stratification in the vertical direction, indicating that the salinity of different layers changes simultaneously with time, with no contribution from salinity convection and diffusion to the bottom layer. A comparison between Figures 6D and F reveals similar vertical distribution characteristics, demonstrating the clear effect of salinity transport from the surface to the bottom layer with time. In Figure 6 , when evaporation is strong, the salinity of the upper surface in panel d is lower than that of panel f, resulting in greater salinity in the deep seawater of the former. This is due to the double effect mentioned earlier, where the strengthening of experimental Ctrl salinity is influenced not only by the sea surface salt flux but also by changes in unit volume after seawater stratification. As a result, the salinity of the deep seawater is enhanced to a certain extent. This effect is also responsible for the results of the NoFlux experiment. However, it is evident that the latter’s impact is smaller than that of the former, and experimental Ctrl can mask the influence of the latter by salt flux to produce an appropriate salinity stratification effect. In contrast, experimental NVDCS only relies on the natural convection and diffusion of salinity for its salinity distribution and transport into the deep sea.

Figures 6A and D clearly show the appearance of a low-salinity stream on the upper surface after 5000 days of simulation in experiment Ctrl. In Figure 9 , it is also evident that the salinity gradient between the high and low latitudes of the sea surface in experiment Ctrl is smaller compared to that in the NVDCS experiment, and the western strengthening effect is more pronounced. As a result, a light saltwater flow moves south from the western boundary, reaching the southern boundary and spreading eastward. Although this phenomenon also exists in NVDCS, its larger salinity gradient makes the uniform western boundary reinforcement effect less effective north of 30°.

At point P2, precipitation and evaporation are balanced, and the vertical distribution at this point is primarily influenced by the exchange with neighboring water bodies. Similar to Figure 6 , Figures 7B and E exhibit no significant stratification. Although experiment Ctrl and NVDCS exhibit comparable structures, a comparison between Figures 7A and B reveals an increase in salinity at the bottom and a decrease in salinity at the sea surface for point P2 in the simulation. The effect is more pronounced in the NVDCS experiment. This is attributed to the transport of low brine from the sea surface in high latitudes towards the south and the northward transport of high brine from the deep layer in low latitudes. The NVDCS model can better simulate this phenomenon compared to other formats.

Point P3 is situated in a high latitude region with intense precipitation. Similar to Figures 6 , 8B, E , 9 do not exhibit any stratification. However, Figures 8D and F clearly demonstrate the vertical salinity stratification state after the addition of fresh water, with Figure 8F indicating lower sea surface salinity compared to Figure 8D . Notably, Figure 8 does not reflect our previous theory, which suggests that experimental Ctrl, located in high latitudes, should have lower deep-sea salinity than experimental NVDCS. The horizontal distribution of seabed salinity after 1, 10, and 20 years of experimental Ctrl and NVDCS operation is depicted in Figure 10 . This figure clearly illustrates that the seafloor salinity is higher in the high evaporation area, while it is lighter in the high precipitation area.

When comparing Figures 9 and 10 , it is observed that both the sea surface and seafloor of Experiment Ctrl had a stronger westward reinforcement compared to Experiment NVDCS. At the surface, Experiment NVDCS had a greater density difference which partly counteracted the north-to-south flow. At the seafloor, Experiment Ctrl responded rapidly to changes in sea surface salinity due to the dual effect of vertical stratification, leading to a greater density difference on the seafloor, which in turn facilitated the flow from south to north. Therefore, Experiment Ctrl demonstrated a stronger westward strengthening effect than NVDCS, regardless of the sea surface or seafloor.

The cross-sectional analysis of Figure 11 shows that both experiments exhibit similar salinity distribution and evolution over time. Freshwater is enriched in high latitudes and subsequently moves southward, leading to a gradual decrease in surface salinity. On the other hand, high saline water is enriched in low latitudes and moves northward via deep transport, leading to an increase in deep seawater salinity and expanding it towards the north. The pressure gradient forces influence the seawater flow, causing low-salt water to move southward. Due to the westward strengthening effect, the flow velocity on the west side increases, eventually leading to the formation of a stream of low-salt water on the southern surface. This observation provides evidence that NVDCS can effectively simulate the vertical salinity distribution without compromising its dynamic process, all while maintaining salinity conservation.

To visually illustrate the differences between the two experiments, we have plotted the disparity of the results from the two experiments in Figure 12 . Positive (negative) values indicate that the salinity in the Ctrl experiment is higher (lower) than in the NVDCS experiment. It is evident that, at the sea surface, the salinity in the Ctrl experiment is higher (lower) in regions of intense precipitation (evaporation) compared to our computational scheme. The opposite trend is observed at the seafloor. This once again confirms the presence of the aforementioned artificial error and demonstrates that our computational scheme effectively addresses this issue.

In the Figure 13 , it can also be seen in the cross-sectional cd that both the Ctrl and NVDCS experiments have a similar structure, in which low-salinity water is concentrated on the east and west sides of the surface layer. The low-salinity water on the western side has a deeper depth, narrower cross-section, and a larger salinity gradient, especially in the Ctrl experiment, which further supports the westward strengthening effect observed in the experiment.

The introduction of freshwater flux, such as precipitation, will alter the mass and salinity of the ocean surface layer, as well as the pressure at the surface and throughout the deep water, thereby affecting the available potential energy. The original calculation method leads to the instantaneous transmission of salinity changes to the seafloor. Meanwhile, the “double effect” results in a greater density difference in the deep water, leading to an overall overestimation of available potential energy, as illustrated in Figure 14 . Equation 34 represents the formula for calculating available potential energy proposed by Benoit and Beckers (2011).