%\affil[2]{Leibniz Institute for Baltic Sea Research Warnem\"nde (IOW), Seestraße 15, 18119 Rostock, Germany}
\begin{document}
\maketitle
\renewcommand{\abstractname}{Abstract}
\begin{abstract}
%
In this article the development and the validation of a high-resolution regional \esm\ is described, to downscale reanalysis or \gcm\ simulations to the Baltic Sea region.
%
Currently, the model consists of the \mom\ for the ocean and the \cclm\ for the atmosphere.
%
This model allows to study regional climate phenomena and to produce climate data appropriate for end users and policy makers.
%
Technically, the model is driven by the atmospheric boundary conditions that are a priory generated from reanalysis data \gcm\ simulation.
%
The data stemming from a run that is forced reanalysis data can be used to validate the quality of the model by comparison to the same reanalysis or observation data.
%
The bi-directional ocean-atmosphere coupling allows allows for a realistic air-sea feedback which clearly outperforms the traditional approach of using uncoupled standalone models as typically pursued with the \EC\ protocol.
%
The coupling of model components is implemented such that is highly scaleable and can potentially be used on \hpc\ centers world-wide.
%
In order to address marine environmental problems (e.g. eutrophication and oxygen depletion), the ocean model encompasses a marine biogeochemistry model setup suitable for the Baltic Sea's hydrographic conditions.
%
The model might be driven by reasonable \ssps\ including different assumptions for nutrient load scenarios.
%
Beside these applications of high societal relevance, the \esm\ can be used for various scientific questions such as climate sensitivity experiments, reconstruction of ocean dynamics, study of past climates and natural variability as well as investigation of ocean-atmosphere interactions.
%
Hence, it can serve for better understanding of natural processes via attribution experiments that relate observed changes to mechanistic causes.
\end{abstract}
\renewcommand{\abstractname}{Kurzfassung}
\begin{abstract}
%
\open{\lipsum[1-2]}
\end{abstract}
\newpage
\tableofcontents
\newpage
\printacronyms
\newpage
\section{Introduction}
\open{\lipsum[1-3]}
\section{Theoretical background and methodology}
One basic problem when dealing with regional climate models system is that the individual components (atmosphere, ocean, land etc.) are described by different models (realized as computer programs) that act on different grids, see \Fig{fig:grids2D} and \Fig{fig:grids1D}.
%
Still, the components have to be coupled in order communicate their state to each other and exchange fluxes as given by nature itself.
Abstraction of figure \ref{fig:grids2D} in one spatial horizontal dimension. Bottom model can support different surface types as water (blue) and ice (white).
}
\end{subfigure}
\caption{The different grids of the different models.}
\end{figure}
\subsection{Introducing the exchange grid and the flux calculator}
\label{sec:exchange_grid_and_flux_calculator}
To exemplify the aforementioned problem let us consider the calculation of a flux of something from the atmosphere to the bottom.
%
Usually, if the flux depends on the state of the ocean, some state variables have to be communicated first to the atmosphere.
%
Since the atmospheric grids are normally larger, the information has to be averaged (weighted by areas) over several bottom cells,
over different surface types, see \Fig{fig:conservative_mapping1}.
With this averaged state information and its own internal state, the atmospheric model can now calculate the flux as a field on the atmospheric grid.
%
It is noteworthy that the flux cannot be calculated for the different surface types differently, since the atmospheric model does not know about that (at least not without changing the code significantly).
%
Finally, the flux field then has to be redistributed on the bottom cells (again in a area-weighted manner such that flux variable is overall conserved),
see \Fig{fig:conservative_mapping2}.
%
Since the flux is only calculated from averaged information and not surface-type-dependent, this approach locally not consistent and can become inaccurate.
%
This is especially true if many bottom grid cells are covered by one atmospheric grid cell.
Calculation of fluxes in the atmospheric model and remapping on to the ocean's grid.
}
\end{subfigure}
\captionsetup{width=\linewidth}
\caption{The standard way of coupling.}
\end{figure}
%
The alternative approach chosen within the developed \esm\ is the introduction of a third component, i.e. the flux calculator that acts on the exchange grid.
%
This grid is the set of intersections between the atmospheric and the bottom grid and thus has, by construction,
a higher resolution than all involved model components, see \Fig{fig:grids1D_plus_exchange}.
Fluxes are calculated and subsequently communicated to the bottom model.
}
\end{subfigure}
\captionsetup{width=\linewidth}
\caption{Coupling the models via the exchange grid and the flux calculator.}
\end{figure}
%
Some fluxes that are not simply determined by surface variables, however, do not fit into this concept.
%
In particular, precipitation and (downward shortwave and longwave) radiative fluxes will still be determined by the atmospheric model and the resulting fluxes will be sent to the flux calculator executable.
%
Subsequently, radiation and precipitation is then simply passed to the bottom models.
%
Still, there is no direct communication between the two physical components and this simplifies ultimately interchangeability of the models.
\subsection{Current implementation}
As the first step, a working version of the coupled \esm\ is developed that consists of the \mom~\citesqr{neumann2021} ocean model and the \cclm~\citesqr{cclm2020} atmospheric model.
%
In contrast to other coupled models, the developed \iow\ \esm\ involves a \textit{flux calculator} executable that mediates the coupling between these models via the so-called \textit{exchange grid}.
%
How the coupling is implemented is described in detail in \Secs{\ref{sec:coupling_cycle} and \ref{sec:flux_formulas}}.
Importantly, as mentioned in \Sec{sec:exchange_grid_and_flux_calculator} the \esm\ is designed such, that other models are added and the current configuration might be extended or replaced by other suitable models.
In the current implementation the \mom\ ocean model and the \cclm\ atmospheric model exchange the following quantities during one coupling time step via the flux calculator, see also \Fig{fig:coupling_cycle}.
%
First, the ocean model sends its state variables, i.e. surface temperature, the albedo and the ice coverage, for each surface type (water and different ice classes; depicted by double arrows in \Fig{fig:coupling_cycle}) and each grid cell to the flux calculator.
%
While the data is transferred it is mapped from the ocean model's grid to the exchange grid.
%
With this input, the flux calculator may compute the \bbr, see \Sec{sec:flux_formulas}, that is emitted by the ocean.
%
This quantity is then sent to both models (and mapped to their grids) where it is added to the atmospheric thermal radiation budget and subtracted from the ocean's one.
%
Note that the thermal radiation that is emitted by the atmosphere is entirely computed in the atmospheric model since it is not simply given as \bbr.
%
Since the atmospheric model also requires the three ocean's state variables mentioned above (for computing transfer coefficients, radiation fluxes and precipitation), they are passed through the flux calculator to the atmospheric model as well.
Second, the atmospheric model calculates its own state variables that are then send to the flux calculator.
%
Using this information the remaining fluxes, i.e. evaporation, latent and sensible heat as well as momentum fluxes, may be computed on the exchange grid and send to both models, see \Sec{sec:flux_formulas}.
%
Note that for the ocean model these fluxes may be calculated depending on the surface type (double arrows in \Fig{fig:coupling_cycle}).
%
Radiation and precipitation fluxes that are not computed by the flux calculator, are simply passed through to the ocean model.
%
Additionally the ocean model requires a few atmospheric state variables, i.e. atmospheric pressure and ten-meter wind-speed components, for \open{...} and the implemented wave model, respectively.
In the current implementation of the \esm, the formulas for the exchanged fluxes are used, that are implemented in the \cclm~\citesqr{cclmManualII2011}.
with the gas constants $R_d$ for dry air and $R_v$ for water vapor.
%
The sea-surface temperature $T_s(x,y,t)$ determines the saturation pressure $p_{\mathrm{sat}}(x,y,t)$ that is calculated according to the Tetens approximation~\citesqr{tetens1930}.
Having the water vapor content $q^s_{v}(x,y,t)$ at hand, one may then calculate the temperature $\tilde{T}$ at which dry air at the surface would show the same energy $p\cdot V$ as the moist air which is there now.
With the density and the coefficients $c_h(x,y,t)$ for turbulent moisture and heat transfer as well as $c_m(x,y,t)$ for the turbulent momentum transfer, one can compute the underlying fluxes according to the following formulas.
%
All fluxes scale with the magnitude of the horizontal wind velocity $\vec{u}(x,y,t)$ present at the water surface.
The evaporation mass flux is calculated assuming that the air adjusts its water vapor content $q_a^v(x,y,t)$ to the one present at the sea surface $q_s^v(x,y,t)$, i.e.
where all quantities are meant to be functions of $(x,y,t)$, i.e. the horizontal location on the sea surface and time, however, for the sake of brevity we skip these arguments.
The \lh\ flux is then directly proportional to the evaporation
%
\begin{align}
\phi_{\mathrm{LH}}(x,y,t) = \Delta H \phi_{\mathrm{evap}},
\end{align}
%
where $\Delta H$ is the constant for the latent heat of either evaporation, freezing or sublimation, depending on the type of phase transition.
The \sh\ flux is determined by the difference between the temperatures of the lowest (discretized) atmospheric layer $T_a(x,y,t)$ and the ocean's surface $T_s(x,y,t)$, i.e.
The momentum fluxes (i.e. the shear stress at the components interface) depend non-linearly on the wind velocity $\vec{u}(x,y,t)$ at the lowest atmospheric layer and are calculated as
where $\xi$ and $u_\xi$ account either for the $x$ or $y$ direction of the velocity field.
%
It is noteworthy, that the horizontal velocity components of the ocean's water body are negligible compared to atmospheric ones.
The thermal radiation that is emitted in upward direction by the ocean is described by the radiation of a black body having the ocean's surface temperature.
%
Thus the thermal flux can be calculated via the Stephan-Boltzmann law suitable for \bbr
%
\begin{align}
\phi_{\mathrm{BBR}}(x,y,t) = \sigma T_s^4,
\end{align}
%
where $\sigma$ is the Stephan-Boltzmann constant.
%
Importantly, since the \bbr\ depends strongly non-linear on the temperature, this flux exemplifies the importance of the local consistency within the coupling.
The downward radiation fluxes (i.e. \sw\ and \lw radiation) as computed by the atmospheric \cclm\ do not only simply depend on surface fields.
%
Thus, these fluxes are still calculated by the atmospheric model and then passed through the flux calculator to the ocean model.
In order to validate the quality of the developed coupled \esm, a setup is used as described in the following section and is illustrated in \Fig{fig:simulation-setup}.
\caption{\label{fig:simulation-setup}\textbf{Model setup for coupled and uncoupled runs.} Atmospheric grid has a resolution of 0.22 by 0.22$^\circ$ whereas the ocean model's grid has the size 3 by 3 \nm.}
\end{figure}
\subsection{Uncoupled atmospheric model}
First, ERA5~\citesqr{era5} reanalysis data is prepared as forcing/boundary data for the \cclm\ atmospheric model for the time range \FullTime.
%
With these forcing files an \text{uncoupled}\cclm\ run is performed over the \EC~\citesqr{jacob2014eurocordex} domain using a resolution of 0.22$^\circ$ by 0.22$^\circ$.
\subsection{Uncoupled ocean model}
From the resulting atmospheric model output, a meteorological forcing for the \mom\ stand-alone ocean model for the Batlic Sea is generated.
%
The uncoupled \mom\ run is then performed with a horizontal resolution of 3$\times$3 \nm.
%
\open{The Baltic Sea's \obcs\ at the Skagerrak stems from ...}
%
\open{The river runoff ... }
%
The marine bio-geochemistry is modeled by the latest version of the internally coupled \ergom~\citesqr{neumann2021}.
\subsection{Coupled \esm}
For the coupled \esm\ run the same atmospheric (ERA5) forcings is used as for the uncoupled \cclm run.
%
Likewise, for the ocean model the same \obcs\ and river runoff data has been used as for the uncoupled \mom\ run.
%
The other input parameters of the model components are also kept the same as in both, atmospheric and ocean model, uncoupled runs.
%
However, the interaction between ocean and atmosphere over the Baltic Sea region is now realized as a two-way coupling via the \textit{exchange grid} and the \textit{flux calculator}.
%
(Thus, the forcings generated from the uncoupled \cclm\ run are not needed for the coupled run.)\\
The data of all three kinds of runs, i.e. uncoupled \cclm, uncoupled \mom\ and the coupled \esm, is analyzed in the following for a time span of 60 years, i.e. \ValidationTime.
%
The employed analysis is described in detail in the following and the resulting figures are shown and discussed in \Secs{\ref{sec:CCLM} and \ref{sec:MOM5}}.
\section{Validation procedure}
In order to judge the quality of the individual model runs the data has been compared to ERA5 reanalysis data.
\caption{\label{fig:flow-diagram}\textbf{Diagram for the data flow during automatized validation procedure.} The main data flow is from left to right and is exemplified for the ocean model's variable \texttt{SST} (sea surface temperature). The arrows depict the particular flow of the data from one validation task to the next. The raw model/reference data (left column) is first pre-processed into a generic format (right column) and subsequently analyzed by \cdo\ operations (upper row). The figures are then plotted by Python tools (lower row) and compiled into a validation report (bottom).}
First, the raw data of one model within the \esm\ as well as the chosen reference data are pre-processed such that all considered model variables are stored in files with the name of the model variable.
For instance the MOM variable \texttt{SST} (sea surface temperature) is stored in file \texttt{SST.nc} that is the broadly used NetCDF~\cite{rew1990netcdf} format.
Practically these steps are implemented as calls of the \cdo~\cite{schulzweida_uwe_2019_2558193} tool and can be configured via a configuration file that is described below.
Note further that the pre-processing might be very specific to the model/reference since usually different models/references feature different data formats.
However, after these very first pre-processing steps the model/reference data format is intended to be model-independent since only the variable's name is needed to perform the following analysis steps.
This generic data format avoids code duplication since all of the analysis and the plotting scripts can then be equally applied to different models within the \esm\ having different reference data.
The analysis and the plotting steps of validation procedure are configured in an own module implemented in the programming language Python~\cite{python3} via a so-called Python \textit{dictionary}.
%
Each model variable of interest is assigned to a key with the variable's name in that dictionary.
%
The value belonging to this key is a dictionary itself, that can contain the
\begin{itemize}
\item seasons for temporal means, i.e. a list months to be averaged over
\item stations defined by coordinates
\item regions defined by coordinates of a rectangle or a mask file that cuts the particular region from the data
\item temporal mean operations that can be applied to time series data and vertical profiles that are extracted for the stations and regions
\item file pattern to the reference data files that contain the corresponding reference variable
The plots are then compiled into a single report that is either provided as a Markdown~\cite{markdown-guide} file, an interactive Jupyter Notebook~\cite{kluyver2016jupyter} or as ready-to-read PDF file.
The Jupyter Notebook format enables the user to adapt the individual plotting scripts to customize figures to be used in publications or other documents.
%
An example of two validation reports in the PDF format can be found in \Sec{sec:appendix}.
In order to investigate the impact of the described exchange grid approach, see \Sec{sec:exchange_grid_and_flux_calculator}, two alternative exchange grid setups are considered for comparison.
%
For the sake of clarity within this discussion, we first have to generalize the exchange grid term as it was introduced above.
In \Sec{sec:exchange_grid_and_flux_calculator}, the term exchange grid was synonymously used for the grid that is formed from the intersections of the involved models grids.
%
In the following we want to refer to this (in fact) \textit{special case} of an exchange grid, as the \textit{intersection grid}.
The two alternative exchange grid cases that are apparent, can then be constructed from either the \textit{atmospheric model grid} or the \textit{ocean model grid} itself.
Since for a typical coupled model setup, the atmospheric grid has the lower resolution than the ocean model, the two resulting alternative exchange grids can differ quite substantially from each other as well as from the intersection grid.
In any case both alternatives will have (by construction) a lower resolution than the intersection-type exchange grid, see again \Sec{sec:exchange_grid_and_flux_calculator}.
First, we employ the approach introduced in this manuscript and calculate fluxes by the flux calculator with state variables locally resolved on the intersection grid and subsequently communicate the fluxes to the models.
%
Second and third, we employ one of the model grids as the exchange grid and calculate fluxes with spatially averaged fields and communicate then the fluxes to the involved models.
These two last cases include also the ``standard'' coupling approach, i.e. using a conservative mapping of state variables from the ocean to the atmospheric model accompanied by the flux calculation via the latter and the communication back to the former~\open{\cite{citation-needed}}.
\caption{\label{fig:remappings-atmos} Atmospheric grid is used as the exchange grid.}
\end{subfigure}
\captionsetup{width=\linewidth}
\caption{\label{fig:remappings}\textbf{Mapping weights for the different kinds of exchange grids.}
%
The atmospheric grid is depicted by the grey lines, the ocean model's grid corresponds to the dark blue lines and the exchange grid is shown with the thin orange lines.
%
Exchange grid cells are additionally filled with transparent orange color.
%
The opaque background colors (corresponding to the color bar green-blue-white) refer to the mean mapping weight contributing to a particular cell on the destination grid.
%
The white numbers show how many cells contribute to that mean.
%
If there is no number given, then only one cell contributes to the particular cell.
%
The columns of each figure correspond to different phases during one coupling time step (from left to right).
%
Left panels depict the sending of state variables from the models to the exchange grid, right panels the communication of fluxes back to the models.
As it can be seen therein, the different mapping matrices for different exchange grids can be distinguished at which phase of the coupling spatial averaging is performed.
For instance in case of the intersection grid, \Fig{fig:remappings-intersection}, the weights are all equal to one when the model's state variables are mapped to the exchange grid, see white grid cell areas in the left panels therein and note the color bar and figure caption.
%
After the fluxes are calculated from the state variables, the fluxes have to be communicated back to the model grids.
%
This mapping naturally involves averaging over several cells and cannot be avoided since the models do eventually feature different grids.
%
This is illustrated by first the background color of the cells accounting for the mean mapping weight and second by the white numbers that count how many cells contribute to the particular destination grid cell.
%
The higher this counter (and, thus, the lower the average mapping weight is) the more information is lost during the conservative mapping from one grid to the other.
%
Importantly, when using the intersection grid as the exchange grid, the state variables are communicated to the flux calculator without any loos of information.
%
Thus, there is no local inconsistency due to any non-linearity of the flux formulas and no errors stemming from the mapping procedure can be amplified by strongly non-linear dependencies of the fluxes.
%
Averaging over several cells only happens when the fluxes are finally communicated to the models.
In contrast, for the other two cases, there is a loss of information, when the state variables are communicated to the exchange grid for the flux calculation, see \Figs{\ref{fig:remappings-bottom} and \ref{fig:remappings-atmos}}.
%
The case when the ocean model's grid is used, see \Fig{fig:remappings-bottom}, seems to be quite similar to the intersection grid case.
%
However, some local information is lost when mapping the atmospheric state variables to the exchange grid before the flux calculation.
%
One can suppose from \Fig{fig:remappings-atmos}, that largest local inconsistencies will occur when the standard approach is employed, i.e. the ocean's state variables are first communicated to the atmospheric grid, fluxes are calculated by the atmospheric model and finally the fluxes are communicated back to the ocean.
When using the atmospheric grid as the exchange grid it was not possible to integrate the coupled model over the whole time period \FullTime.
%
Instead the model becomes instable after 12 months featuring an unrealistic low \sst\ at specific points on the ocean's grid, see \Fig{fig:SST-crash}, where a snapshot directly before the model stops is depicted, i.e. January 9, 1961.
%
For instance the cell located at 21.08$^\circ$ east and 59.2$^\circ$ north exhibits a temperature that falls rapidly down by 40$^\circ$C within approximately 6 hours, see \Fig{fig:time-series-crash} blue curve.
%
Before this event occurred, the evaporation increases, see orange curve therein.
%
Due to the low exchange grid (given by the atmospheric model, orange boxes in \Figs{\ref{fig:SST-crash} and \ref{fig:EVAP-crash}}) the magnitude of the evaporation is mainly given by the surrounding ice-free cells, see white boxes in \Fig{fig:SST-crash}.
%
Thus, the ice covered cell is cooled down by a rate that is determined by the liquid water contained in the ice-free cells.
%
Hence, the instability is a direct consequence of the inconsistency when calculating fluxes on the low-resolution atmospheric grid.
\subsection{Intersection grid vs. ocean model grid}
\caption{\label{fig:T_2M_AV_anomalies}\textbf{Seasonal anomalies for the two-meter air temperature.} The anomaly is with respect to the ERA5 reanalysis data. The temporal averages are performed over a time range \ValidationTime.}
\end{figure*}
\paragraph{Time series}
See \Fig{fig:T_2M_AV_time_series_stations} and \Fig{fig:T_2M_AV_time_series_regions}.
\caption{\label{fig:T_2M_AV_time_series_stations}\textbf{Time series for the two-meter air temperature at chosen stations.} The anomaly is with respect to the ERA5 reanalysis data. The temporal averages are performed over a time range \ValidationTime.}
\caption{\label{fig:T_2M_AV_time_series_regions}\textbf{Time series for the two-meter air temperature at chosen regions.} The anomaly is with respect to the ERA5 reanalysis data. The temporal averages are performed over a time range \ValidationTime.}
\end{figure*}
\paragraph{Taylor diagrams}
See \Fig{fig:T_2M_AV_taylor_diagrams_stations} and \Fig{fig:T_2M_AV_taylor_diagrams_regions}.
\caption{\label{fig:T_2M_AV_taylor_diagrams_stations}\textbf{Taylor diagrams for the two-meter air temperature at chosen stations.} The anomaly is with respect to the ERA5 reanalysis data. The temporal averages are performed over a time range \ValidationTime.}
\caption{\label{fig:T_2M_AV_taylor_diagrams_regions}\textbf{Taylor diagrams for the two-meter air temperature at chosen regions.} The anomaly is with respect to the ERA5 reanalysis data. The temporal averages are performed over a time range \ValidationTime.}
\end{figure*}
\paragraph{Cost functions}
See \Fig{fig:T_2M_AV_cost_functions_stations} and \Fig{fig:T_2M_AV_cost_functions_regions}.
\caption{\label{fig:T_2M_AV_cost_functions_stations}\textbf{Cost functions for the two-meter air temperature at chosen stations.} The anomaly is with respect to the ERA5 reanalysis data. The temporal averages are performed over a time range \ValidationTime.}
\caption{\label{fig:T_2M_AV_cost_functions_regions}\textbf{Cost functions for the two-meter air temperature at chosen regions.} The anomaly is with respect to the ERA5 reanalysis data. The temporal averages are performed over a time range \ValidationTime.}
\caption{\label{fig:SPEED_10M_AV_anomalies}\textbf{Seasonal anomalies for the ten-meter wind speed.} The anomaly is with respect to the ERA5 reanalysis data. The temporal averages are performed over a time range \ValidationTime.}
\caption{\label{fig:DAY_PREC_anomalies}\textbf{Seasonal anomalies for the precipitation.} The anomaly is with respect to the ERA5 reanalysis data. The temporal averages are performed over a time range \ValidationTime.}
\end{figure*}
\open{
\paragraph{Further plan}
\begin{itemize}
\item now only results from coupled model are shown, next step is to compare to uncoupled (data is available)
\caption{\label{fig:SST_anomalies}\textbf{Seasonal anomalies for the sea-surface temperature.} The anomaly is with respect to the ERA5 reanalysis data. The temporal averages are performed over a time range \ValidationTime.}
\caption{\label{fig:temp_profiles_stations}\textbf{Vertical temperature profiles at chosen stations in the Baltic Sea.} The reference is the BED data set\open{~\citesqr{???}}. The temporal averages are performed over a time range \ValidationTime.}
\caption{\label{fig:FI_anomalies}\textbf{Seasonal anomalies for the sea-ice coverage.} The anomaly is with respect to the ERA5 reanalysis data. The temporal averages are performed over a time range \ValidationTime.}
\end{figure*}
\open{
\paragraph{Further plan}
\begin{itemize}
\item now only results from coupled model are shown, next step is to compare to uncoupled (data is available)
\subsection{Simple model for the cold-bias in the coupled model}
In order to explain the strong cold bias during summer time in the coupled model simulation, a simplistic model is considered.
\subsubsection{The simplistic model}
The atmospheric temperature directly above the water surface $T_a$ and the water-surface temperature $T_w$ are written as functions of a set of quantities, i.e.
where $\{A_i\}$ and $\{W_i\}$ are all \dofs\ in the atmosphere and the water, respectively, that are not of interest now.
%
The downward shortwave radiation flux that heats up the ocean is explicitly taken into account via $\phi_{\sw}$ as well as the latent heat flux $\phi_{\lh}$ that mainly cools the ocean.
%
The total change of these quantities can then be written down as
The different signs are due to the assumption $\gamma > 0$ and $\delta < 0$.
%
Consequently, in the coupled case the impact of a changing shortwave radiation and latent heat flux on the water's temperature is \textit{amplified} with respect to the uncoupled case.
%
Thus, reducing the incoming shortwave radiation ($\d\phi_{\sw} < 0$) and increasing the latent heat flux ($\d\phi_{\lh} > 0$) we decrease the water's temperature more strongly in the coupled case.
%
Hence, a deficit in the radiation flux and an overestimated latent heat (as it is present for CCLM) leads to a stronger cold bias for the coupled model.
