A Moist Quasi-Geostrophic Coupled Model: MQ-GCM2.0

A Moist Quasi-Geostrophic Coupled Model: MQ-GCM2.0 Sergey Kravtsov1,2,3, Ilijana Mastilovic1, Andrew McC. Hogg4, William K. Dewar5,6 and Jeffrey R. Blundell7 1Department of Mathematical Sciences, University of Wisconsin-Milwaukee, P. O. Box 413, Milwaukee, WI 53201, USA 2Shirshov Institute of Oceanology, Russian Academy of Sciences, Moscow, 117218, Russia 5 3Institute of Applied Physics, Russian Academy of Sciences, Nizhniy Novgorod, 603155, Russia 4Research School of Earth Sciences, and ARC Centre of Excellence in Climate Extremes, Australian National University, Canberra, Australia. 5Department of Earth, Ocean and Atmospheric Science, Florida State University, Tallahassee, FL 32304, USA 6Laboratoire de Glaciologie et Geophysique de l'Environnement, CNRS, Grenoble, France 10 7Ocean and Earth Science, National Oceanography Centre Southampton, University of Southampton, Southampton SO14 3ZH, United Kingdom

documented in the Q-GCM Users' Guide, v1.5.0 (Hogg et al., 2014). The model couples the multi-layer quasi-geostrophic 35 (QG) ocean and atmosphere components via ageostrophic mixed layers that regulate the exchange of heat and momentum between the two fluids. Q-GCM model can be configured as either a box (a 'double-gyre') or a channel ocean (a 'Southern Ocean') underneath a channel atmosphere; it conceptualizes the mid-latitude climate system driven by the latitudinal variation of the incoming solar radiation. In addition to the oceanic mixed layer, the model physics incorporates a dynamically active atmospheric mixed layer (effectively, the atmospheric planetary boundary layer: APBL), the dependence 40 of the wind stress on the ocean-atmosphere surface velocity difference, as well as a dynamically consistent parameterization of the entrainment heat fluxes between the model layers. It can also be easily modified to include a parameterization of a seasurface temperature (SST) feedback on the wind stress (e.g., Hogg et al., 2009), which will be a part of the new version of the model developed here. Q-GCM thus encompasses a richer, more comprehensive set of processes, enabling one to achieve a more accurate simulation of the ocean-atmosphere coupling, especially at mesoscales, relative to some other, analogous 45 conceptual models, which either assume the atmospheric near-surface temperature to be in equilibrium with SST (e.g., Feliks et al., 2004cf. Schneider and Qiu, 2015) or relate this temperature in an ad hoc way to the instantaneous in-situ distribution of the model's tropospheric temperature (Kravtsov et al. 2005(Kravtsov et al. , 2006Deremble et al. 2012). The Q-GCM model was previously used for ocean-only and coupled experiments around the double-gyre problem (Hogg et al., 2005Martin et al. 2020), as well as in the ocean-only studies of the Southern Ocean's climate system (Hogg and 50 Blundell, 2006;Meredith and Hogg, 2006;Hogg et al., 2008;Hutchinson et al., 2010;Kravtsov et al., 2011).
The long oceanic thermal and dynamical inertia makes the ocean a primary agent for generating potentially predictable climate signals on time scales from years to decades, whereas atmospheric intrinsic time scales are significantly shorter. The null hypothesis for climate variability views the ocean as a passive integrator of high-frequency noise associated with atmospheric geostrophic turbulence (Hasselmann, 1976;Frankignoul and Hasselmann, 1977;Frankignoul, 1985;Barsugli 55 and Battisti, 1998;Xie, 2004). However, both observations (e.g., Chelton, 2013;Frenger et al., 2013) and decades of experimentation with wind-driven eddy-resolving ocean models (e.g., Berloff and McWilliams, 1999;Primeau, 2002;Hogg et al., 2009;Shevchenko et al., 2016, among many others) documented vigorous internal variability and the associated mesoscale features (fronts and eddies with spatial scales of 10-100km) throughout the World Ocean. Two key regions in which these eddies are most important are near the western boundary currents and their 60 extensions and in the Southern Ocean. Mesoscale variability in these regions modulates atmospheric fronts and storms' intensity and distribution, thus affecting atmospheric variability on short time scales (e.g., Maloney and Chelton, 2006;Minobe et al., 2008;Nakamura and Yamane, 2009;Bryan et al., 2010;Chelton and Xie, 2010;Kuwano-Yoshida et al., 2010;O'Neill et al., 2010O'Neill et al., , 2012Frenger et al., 2013;O'Reilly and Czaja, 2015;Seo et al., 2016;Parfitt et al., 2017). Recent observational and modelling evidence strongly suggested that this mesoscale oceanic turbulence may also imprint itself onto 65 large-scale low-frequency climate modes (with time scales from intra-seasonal to decadal), which would have profound consequences for near-term climate predictability (e.g., Siqueira and Kirtman, 2016). To study this https://doi.org/10.5194/gmd-2021-160 Preprint. Discussion started: 18 June 2021 c Author(s) 2021. CC BY 4.0 License. phenomenon may, therefore, require coupled climate models with high horizontal resolution in both their oceanic and atmospheric components (see below); this would make the requisite long climate simulations using highly resolved state-ofthe-art climate models computationally infeasible. 70 Feliks et al. (2004Feliks et al. ( , 2007Feliks et al. ( , 2011 and Brachet et al. (2012) examined the response of the atmosphere to mesoscale sea-surface temperature (SST) anomalies through hydrostatic pressure adjustment in an idealized atmospheric model. They showed that resolving an ocean front and mesoscale eddies affects atmospheric climatology, intraseasonal modes, as well as decadal variability (when forced with the observed SST history) in their model (see also Nakamura et al., 2008). These authors argued that atmospheric components of global climate models must resolve oceanic fronts to faithfully simulate the observed 75 climate variability (see also Minobe et al., 2008). Ma et al. (2017) also concluded that "It is only when the (atmospheric) model has sufficient resolution to resolve small-scale diabatic heating that the full effect of mesoscale SST forcing on the storm track can be correctly simulated," with the ensuing consequences for atmospheric low-frequency variability associated with the downstream Rossby wave breaking (Piazza et al., 2016) and blocking (O' Reilly et al., 2015). By contrast, Bryan et al. (2010) proposed that accurate representation of mesoscale ocean-atmosphere coupling in a model depends more on the 80 Marine Atmospheric Boundary Layer (MABL) mixing scheme than on the ability of an atmospheric model to resolve a thermal front per se; this would justify the use of atmospheric resolutions on the order of 50 km in many GCM studies that documented a pronounced influence of the mesoscale air-sea interactions on the atmospheric storm tracks (Miller and Schneider, 2000;Nakamura et al., 2008;Taguchi et al., 2009;Kelly et al., 2010;Perlin et al., 2014;Small et al., 2014;Ma et al., 2015Ma et al., , 2017Piazza et al., 2016). 85 A numerically efficient intermediate-complexity Q-GCM model thus provides an alternative (to highly resolved GCMs) and unique tool ideally suited to help advance our understanding of multi-scale ocean-atmosphere interactions and their climatic impacts. Its QG dynamical core resolves well the geostrophic turbulence on either side of the ocean-atmosphere interface, including oceanic mesoscale eddies/fronts and atmospheric storm tracks. The existing version of the coupled model, however, lacks the parameterization of SST effects on the model's MABL winds. One such parameterization was tested in 90 the Q-GCM's ocean-only configuration by Hogg et al. (2009). Mastilovic and Kravtsov (2019) examined the effects of both Hogg et al. (2009) and Feliks et al.'s (2004 SST-dependent MABL wind formulations in the context of coupled Q-GCM simulations with the standard (coarse) and fine (mesoscale resolving) atmospheric grid spacing. They found thatconsistent with previous studies (Dewar and Flierl, 1987;Maloney and Chelton, 2006;Hogg et al. 2009;Gaube et al., 2013Gaube et al., , 2015Chelton, 2013;Small et al., 2014) -these effects constitute a negative feedback on the ocean and tend to reduce the 95 intensity of the oceanic mesoscale perturbations that generated the mesoscale wind anomalies in the first place. In a coupled setting, this leads to a dilution of oceanic mesoscale features and the resulting lack of model's sensitivity to atmospheric resolution. Surprisingly, the atmosphere-only Q-GCM simulations with and without an SST front, or with and without SSTdependent AMBL winds, apparently also produce statistically identical atmospheric variability irrespective of the atmospheric resolution (Mastilovic and Kravtsov 2020, personal communication), in sharp contrast to the resolution-100 dependent dynamics documented in Feliks et al. (2004Feliks et al. ( , 2007. Furthermore, the entire previous experience with Q-GCM indicates the absence, in the existing version of the model, of a nonlinear, "weather-regime"-type atmospheric behaviour documented in analogous atmospheric (Marshall and Molteni, 1993;Kravtsov et al., 2005) and coupled models (Kravtsov et al., 2006; such behaviour may lead to a nonlinear atmospheric sensitivity to ocean induced SST anomalies and generate fundamentally coupled decadal climate modes (e.g., Kravtsov et al., 2008). 105 The purpose of this paper is to document a new, revamped version of the Q-GCM model, which addresses an apparent lack of the Q-GCM atmosphere's sensitivity to SST anomalies, contrary to ample theoretical and observational evidence otherwise. In section 2, we briefly summarize the dynamical (QG) core of the model, placing some of the supporting information in the appendix. Major modifications to the original, dry-model physics (section 3) include an improved radiative-convective scheme resulting in a more realistic model mean state and the associated model parameters, a new 110 formulation of entrainment in the atmosphere, which prompts a more efficient communication between the atmospheric mixed layer and free troposphere, as well as an addition of temperature-dependent wind component in the atmospheric mixed layer and the resulting mesoscale feedbacks. The most drastic change is, however, the inclusion, in the model, of moist dynamics (section 4), which may be key to midlatitude ocean-atmosphere coupling (Czaja and Blunt, 2011;Laîné et al., 2011;Deremble et al., 2012;Willison et al., 2013;Foussard et al., 2019). Accordingly, this version of the model is to be 115 referred to as the MQ-GCM model. Overall, the MQ-GCM model is shown to exhibit a rich spectrum of behaviours reminiscent of many of the observed properties of the Earth's climate system (section 5). The paper concludes with some discussion in section 6 and section 7, which summarizes the MQ-GCM changes and code modifications with respect to the original Q-GCM version. This presentation is also supplemented with the collection of Fortran 90 routines containing all the new MQ-GCM source code that complement the original Q-GCM distribution (http://www.q-gcm.org). 120 2 Q-GCM dynamical core Q-GCM model incorporates quasi-geostrophic dynamics on a -plane in its n-layer oceanic and atmospheric modules (below, we will use n=3); these dynamics are governed by the equations describing the evolution of quasi-geostrophic potential vorticity . For a flat-bottom ocean (used here for simplicity, although the topography is included in Q-GCM) 125 ( 2) and (3) Note that the layer indexing in the equations above goes from the surface upward and the Laplacian friction term is omitted; otherwise, these equations are completely analogous to the oceanic equations (the atmospheric variables are denoted above by the left superscript 'a'). At the lower atmospheric boundary, the entrainment flux solely represents, in the 150 original Q-GCM formulation, Ekman dissipation (thus signifying the momentum transfer from the atmosphere to the ocean); in the modified version of the model to be developed here it will also include a temperature dependent component capable of driving mesoscale air-sea interaction (section 3.3). The atmospheric model (and thus the entire coupled model) is driven through interior entrainment fluxes in Eq. (4) ( is set to zero), which are a by-product of perturbing the meanstate radiative-convective equilibrium (section 3.1) by a latitudinally non-uniform insolation. The new radiation/heat 155 exchange formulation and mixed layer/entrainment formulation developed here (section 3) are aimed to help achieve a parameter regime with enhanced (and, arguably, more realistic) coupling between oceanic and atmospheric dynamics in the model. They are further modified in formulating the new version of the model with an active hydrological cycle and the associated latent-heat feedbacks (section 4).

Updates to the original, 'dry' version of Q-GCM
In describing the updates below, we will generally focus on the elements of Q-GCM model that have been revised here and quote the values of new parameters, or the updated values of the previously used parameters along the way, while referring the reader/user to the existing Q-GCM guide (Hogg et al., 2014) for a more thorough description of the default model configuration. 165

Radiative-convective equilibrium, atmospheric mean state and convective fluxes
The previous version of Q-GCM assumes purely radiative equilibrium to compute the atmospheric mean state. In the revised version, this assumption is replaced by that of the radiative-convective mean-state balance. We denote the actual (not potential) vertically averaged temperatures within each of the interior atmospheric layers as , 1, 2, 3 to write, over ocean, 170 Here the dot denotes the time derivative and other notations follow the Q-GCM Users' Guide, v1.5.0 (Hogg et al. 2014). In particular, the upward/downward arrow subscripts denote upwelling/downwelling longwave radiative fluxes within each of a e 0 = a w ek a e 1 , a e 2 a e 3 the layers k = m (mixed layer), 1, 2, 3 or the surface (0), the subscripts e-and e+ denote entrainment fluxes below and above interface k, the subscript s refers to the solar radiation and -to the ocean-atmosphere sensible/latent heat exchange. The 175 fluxes describing non-radiative heat exchange between the layers are interpreted, in the mean state, as convective fluxes parameterized in the following way (cf. Manabe and Strickler 1964;Manabe and Wetherald 1967;Ramanathan and Coakley 1978): The upward fluxes in Eqs. (8) are all positive (or are otherwise set to zero), with the coefficient W m -2 K -1 and the 180 critical lapse rate K/km. The potential temperatures of the atmospheric layers are given by (9) where is the dry adiabatic lapse rate (of about 10 K/km).
The radiative fluxes in Eqs. (7) are parameterized assuming that the atmospheric layers have constant emissivity , and the Stefan-Boltzmann expressions for perturbation fluxes are linearized with respect to the basic-state: 185 where (11) 1, 2, 3 and is the Stefan-Boltzmann constant.
To solve for the mean state, we set all of , , , , , to zero and numerically integrate Eqs. (7)-(11) to 190 equilibrium, using Euler differences in time with the time step of 5 min. Setting along with W m -2 results in the mean state whose parameters are listed in Table 1. The atmospheric optical depth decreases with altitude, but so do the unperturbed thicknesses of our chosen atmospheric layers, making the constant layer emissivities above a reasonable first approximation commensurate with an idealized nature of the present model. The model has a realistic (time-mean, global-mean) vertical temperature distribution. Note that the atmospheric reduced gravities are derived, 195 in this version of the model, from the mean-state parameters rather than being prescribed (at 1.2 and 0.4 ms -2 ), as in the previous version (see the appendix for further details). The climatological solution above is formally obtained over ocean, but it also applies over land of zero heat capacity (since the steady state does not depend on the heat capacity of the surface); the land's zero heat capacity is also a feature of the original Q-GCM formulation. The near-surface convective fluxes, however, would generally be different over ocean and over land (which occupies a significant fraction of the atmospheric 200 channel, including fairly large strips both north and south of the ocean to avoid the distortion of ocean-atmosphere interaction by the effects associated with the atmospheric boundary conditions); the values of the convective fluxes in Table   F  1 should thus be interpreted to represent zonally averaged fluxes. Below, in section 3.2, we will describe, among other things, modifications of the atmospheric mixed layer (a.m.l.) perturbation equation (that is, the one describing evolution of the anomalies with respect to the mean state) over land regions. 205

Mixed-layer perturbation equations and entrainment formulation
The perturbation equations in the mixed layers, for primed variables, have the same form as the first two equations (7), aside from addition of advective and entrainment fluxes which we will discuss further below; hereafter, we will drop primes in all perturbation equations for convenience. We will assume that the atmospheric perturbation temperature is vertically uniform, 210 consistent with an active role of convection processes (cf. Manabe and Wetherald, 1967), viz.
(12) (∆ ! " above denotes the potential temperature jump across the k-th interface; see Table 1) which allows one to express all perturbation radiative fluxes in (10) via perturbation oceanic and atmospheric mixed layer temperatures and interfacial displacements ; in particular, 215 (13) where all of the A, D, E coefficients can be written in terms of the known mean-state parameters. The equations (13) above should be compared with (4.2-4.6) of the Q-GCM User's Guide (Hogg et al. 2014). Notably, the parameter was effectively set to 1 in the previous version of Q-GCM, and hence all of the E coefficients were equal to zero. On the other hand, that previous version had additional parameters B and C for radiative corrections associated with the variable a.m.l. 220 depth and topography. Here we are back to the model with a constant a.m.l. depth (see below); we also neglect topography corrections for simplicity.
Another consequence of the assumption , used here, is that the coefficients A, D, E in Eqs. (13) over ocean and over land are different, and so is the a.m.l. temperature equation -the second equation in (7). In particular, over land we have (neglecting, for now, advection and entrainment terms in the a.m.l. equation) 225 (14) where is the infrared upward flux from the surface of the land, the latter assumed to have zero heat capacity [hence zero on the left-hand side of the first equation in (14)] and conductivity (hence ). From Eqs. (14) it follows that over land (15) while, in analogy with Eqs. (10) 230 (16) [compare this with the second Eq. (13)]. The first (additional) term in Eq. (16) will also modify all other upwelling radiation fluxes ( ) accordingly through Eqs. (10), resulting in modified values of the A and D coefficients and all E coefficients set to zero over land.
Yet another, minor, but potentially fairly important modification of the previous Q-GCM formulation is the inclusion of the 235 dependence on the relative wind speed in the bulk formulas for the sensible/latent heat ocean-atmosphere exchange, which plays a significant role in setting up the North Atlantic SST tripole variability (Deser and Blackmon, 1993;Kushnir, 1994;Czaja and Marshall, 2001;Kravtsov et al., 2007;Fan and Schneider, 2012): cf. Eqs. (4.30-31) of Hogg et al. (2014). The entrainment in the ocean interior only occurs between layers 1 and 2, and is computed in the same way as in the original model:  [Hogg et al. 2014, Eq. (4.32)]; following, again, Hogg et al. (2014), is also corrected to have zero area integral by adding a spatially uniform offset value at each time step.
Similarly, in the atmospheric mixed layer (21) Here is the perturbation temperature given by Eq. (12); in the mean state this temperature is set to . Using 260 instead of in Eqs. (21) is what keeps the instability described by Hogg et al. (2003) in check in the present version of the model with constant . This is due to the fact that is tied to the instantaneous vertical structure of the atmosphere, which limits the magnitude of entrainment heat fluxes (as tends to be closer to than ) and also provides additional Hogg et al. (2014) assumed We modify this assumption by making the entrainments across the lower and upper atmospheric interface be linearly related, with the coefficient f2 (see below); this procedure can also be adapted for the use in an n-layer model by introducing additional free parameters analogous to f2. To allow one a degree of freedom in controlling damping rates at each interface somewhat independently, we also introduce here a (small) vertical diffusion, using a 280 linearized version of McDougall and Dewar (1998) formulation: We can now solve the system (23)-(24) for the two unknown non-diffusive entrainment rates and and, hence, for the full entrainment rates and . We use the nominal value of 0.0001 m s -1 for both of and and initially set (25)  285 to ensure generation of similar velocity shears between atmospheric layers 1/2 and 2/3 by the thermal forcing of a given amplitude. Increasing would tend to increase the geostrophic zonal velocity shear between the lower two atmospheric layers and decrease the velocity shear between the upper two atmospheric layers; setting 0 recovers the previous Q-GCM formulation. The optimal value of is to be determined by trial-and-error tuning of the model.

Temperature-dependent flow in atmospheric mixed layer and partially coupled setup 290
We introduce temperature dependence of the a.m.l. winds by modifying the mixed-layer momentum equations in two ways, namely: (i) including, explicitly, temperature driven pressure gradients (which takes into account the mixed-layer hydrostatic adjustment to temperature contrasts: Lindzen and Nigam, 1987), following Feliks et al. 2004Feliks et al. , 2007Feliks et al. , 2011and (ii) making the surface drag coefficient depend on the ocean-atmosphere temperature difference to parameterize changes in a.m.l. stability (Wallace et al., 1989), following Hogg et al. (2009); see Small et al. (2008) and Chelton and Xie (2010) for a review 295 of these two mechanisms for mesoscale air-sea coupling. Putrasahan et al. (2013) demonstrated that, in the Kuroshio region, both mechanisms (i) and (ii) are important, with relative contributions depending on the spatial scale of the SST anomalies. Putrasahan et al. (2017) also concluded that heat advection by oceanic mesoscale currents plays a key role in creating such SST anomalies and forcing the MABL response in the Gulf of Mexico. To implement these changes, we write the a.m.l. momentum equations as 300 (26) where is the geostrophic velocity in the lowest atmospheric layer, is the wind stress, , =1 and =0.15. Setting one of the -parameters to zero can be used to examine processes (i) and (ii) above independently; setting both of these parameters to zero would recover the previous, temperature-independent a.m.l. wind formulation (3.2)-(3.3) of Hogg et al. (2014). On top of these modifications, we also set the drag coefficient over ocean to 2/3 of the default 305 value over land, following Marshall and Molteni (1993).
Upon adding to (26) analogous equations for oceanic mixed layer in their original form [Hogg et al. 2014;Eqs. (3.4)-(3.5)], we end up with a closed system of equations for the unknown values of at each grid point, which can be solved analytically in the same way as before (see Hogg et al., 2014). Note that additional temperature gradients in the first two equations (26) produce a non-divergent wind field with zero direct Ekman pumping and, also, zero temperature 310 advection contributions; their dynamical effect is thus purely indirect, via modifications to the wind-stress field; they also generate non-zero moisture advection in the moist version of the Q-GCM model, to be developed later in section 4.
The temperature-dependent a.m.l. wind formulation (26) is associated with coupled feedbacks that tend to suppress oceanic turbulence and SST fronts (cf. Hogg et al., 2014; see also section 5). In principle, realistic mesoscale ocean field can still be achieved in inherently more turbulent oceanic regimes at high Reynolds numbers, but this requires very high ocean 315 resolution and is computationally demanding. An alternative fix is to apply partial momentum coupling of the oceanic and atmospheric mixed layers in which the atmosphere "sees" the wind stress as per the full version of Eqs. (26), whereas the oceanic wind stress is computed from Eqs. (26) in which In this way the mesoscale feedbacks of temperaturedependent wind, which damp oceanic turbulence, are artificially suppressed, but their effect on the atmosphere is preserved, possibly leading to coupled dynamics involving large-scale low-frequency reorganization of the wind field and the ensuing 320 ocean response.

Lateral boundary conditions for mass and temperature equations
The original Q-GCM formulation employed no-through-flow conditions on the zonal boundaries of the atmospheric channel but effectively allowed the mass to leave/enter ocean mixed layer through side boundaries to avoid Ekman-pumping 325 singularities there [via the direct use of Eqs. .

Moist version of the model: MQ-GCM 340
Perhaps the most important change to the original Q-GCM formulation is the inclusion of the hydrological cycle and latentheat feedbacks, resulting in what we refer to as the Moist Quasi-Geostrophic Coupled Model (MQ-GCM). Indeed, Czaja and Blunt (2011) proposed that the oceans can influence the troposphere through moist convection over the regions with strong mesoscale variability; see also Willison et al. (2013). To compute moisture variables in the model, we assume that the vertical temperature profile at a given (x, y) location is linear in z, temperature decreasing with altitude z above the sea level 345 at the critical lapse rate : Here is the absolute temperature (in K) in layer k [k can be a symbol (m) when referring to the a.m.l. temperature or an index (k=1, 2, 3) when denoting the interior (quasi-geostrophic) layers of the atmospheric model], while in the interior are given by Eq. (12). In such a constant-lapse-rate atmosphere, the pressure p(z) is related to temperature as 350 (29) where g is the gravity acceleration and R is the ideal gas constant for dry air. Combining the latter two equations and the ideal gas law at the basic state with , we can compute the representative densities of each layer by estimating them at the altitude z corresponding to the mid-layer height (e.g., for the mixed layer, for layer 1, etc.); this gives, for the parameters in Table 1 and hPa, 355 kg m -3 . The saturation specific humidity is given by (30) where is the ratio of the dry-air and water-vapor gas constants and the saturation water vapor pressure is computed as in Bolton (1980), using the parameters hPa, , , K and K.
Given the a.m.l. perturbation temperature and the atmospheric interface displacements , the equations (12), (28) Here are the geostrophic velocities in the atmospheric layer k, E is the evaporation and is the precipitation; and are the moisture entrainment fluxes below and above the interface k, respectively (all of these fluxes are in kg 370 m -2 s -1 ). Once again, k= (m, 1, 2, 3) and for . The equations (31) also use boundary conditions analogous to those for temperature (section 3.4).
The evaporation over the ocean is given by (Gill 1982) where the atmospheric specific humidity near the ocean surface is computed assuming constant relative humidity in the 375 a.m.l., following Deremble et al. (2012); the coefficient . Over land, we specify the (fixed in time) evapotranspiration flux (which also includes the zonally averaged evaporation from other ocean basins absent in our onebasin configuration; this allows us to achieve reasonable values characterizing the moist model's climatological distribution of specific humidity). In space, this y-dependent flux decreases linearly from the maximum value of 1 m year -1 at the southern boundary of the atmospheric model (note the usage of water density =1000 kg m -3 here, leading to precipitation 380 estimates in terms of the equivalent water depth per unit time) to the minimum value of 0.1 m year -1 at the northern boundary.
Entrainment fluxes of moisture in Eqs. (31) are formulated in a way analogous to the entrainment heat fluxes. In particular, at the top of the a.m.l. forward with all the precipitation rates set to zero to update the values of the specific humidity; recall that the specific humidity is assumed to be independent of z in each layer. Then the vertical integrals of the newly computed specific 395 humidity excess over the saturated specific humidity (which is a function of z) within each layer are computed [for numerical efficiency, this is done semi-analytically by fitting a quadratic function of z to hk-hs(z)]. This amount of moisture is set to fall out, over the period associated with the leap-frog time step, as the precipitation , and the specific humidity in the corresponding layer is reduced accordingly.
The hydrological cycle above is coupled with the model dynamics via the associated latent heat exchange/release. In the 400 MQ-GCM, the equation (17) is only meant to describe the sensible heat exchange between the ocean and the atmosphere, with a reduced value of the sensible-heat exchange coefficient . In addition, the oceanic mixed layer is experiencing the (perturbation) latent heat loss (in W m -2 ) of (35) and the atmospheric layers -the (perturbation) latent heat gain of 405 where J kg -1 is the latent heat of vaporization of water and k= (m, 1, 2, 3). Note that the full latent-heat fluxes both include the part associated with the basic state of the model in its radiative-convective balance, but the Q-GCM is formulated as a perturbation model, which requires the subtraction of the basic-state latent-heat fluxes. We here assume that the basic-state part of and is approximately given by the spatial averages of over the oceanic 410 basin and atmospheric channel, respectively (this assumption is justified post hoc by the moist model's a.m.l. and o.m.l. climatological temperatures being close to those of a dry model) -and ; hence, we remove these spatial averages at each time step in Eqs. (35) and (36) to define the latent-heat flux anomalies that force our perturbation heat equations.
The fluxes and directly enter the right-hand side of the o.m.l. and a.m.l. equations (18), respectively. The interior latent-heat release is added to the right-hand side of the corresponding layer's heat equation in Eqs. (22), so 415 that the sum enters, as an additional term, the right-hand side of Eq. (23) and modifies the entrainment rates .
Note again, however, that the moisture equations (31) use the first-guess, unmodified ("dry") entrainment rates computed from the original Eqs. (23), (24), upon which the precipitation rates and latent heat corrections are computed and used to adjust the entrainment rates as follows (37)  420 These modified entrainment rates are then used to timestep the atmospheric QGPV equations.
The above numerical scheme -with the first guess 'dry' entrainment driving the moisture equations to produce latent-heatdriven corrections to the entrainment, which are then used to force the QG model interior, -can be further improved by iterating the solution of the moisture equations at a given time step to achieve mutually consistent estimates of both precipitation and entrainment in the interior QG layers. In this scheme, the interior entrainment fluxes at a given iteration 425 would be used, along with the fixed advective and diffusive fluxes, to update the interior humidity and compute the precipitation rates until these rates (and entrainment rates) converge to a steady solution. This procedure is, however, much more numerically challenging than its first guess 'dry entrainment' implementation. The latter 'dry entrainment' ∑ a e k P k a Q L,k a e 1 → a e 1 + Δ a e 1 ; e 2 → a e 2 + f 2 Δ a e 1 ; implementation may formally be justified if the "moist" corrections to the dry entrainment produce relatively small changes to the interior precipitation. To explore this issue further, we included, in the present version of MQ-GCM model, an option 430 which allows one to modify the 'dry' entrainment based precipitation estimates via one additional iteration in which the interior moisture equations (31) are stepped with the entrainment moisture fluxes (34) utilizing the moisture corrected entrainment rates (37).
Moisture transport and latent heat release driving the atmospheric response in the areas away from the oceanic warm regions of evaporation are important elements of air-sea interaction over variable SST fronts, which are altogether missing in a dry 435 version of the Q-GCM model (cf. Deremble et al. 2012).

Model simulations
We set the amplitude of the variable incoming radiation to 150 W m -2 (as compared to 80 W m -2 in Hogg et al., 2014) and ran three 130-yr simulations of the new dry version of the model, as well as of the new moist model, MQ-GCM (six simulations total), using model setups with and without temperature dependence in the atmospheric mixed-layer wind. 440 For both dry and moist versions of the model, the control run, without this dependence, was started from the final state of the preliminary 100-yr spin-up simulation (from rest), the partially coupled simulation was initialized by the final state of the control run, and the fully coupled temperature-dependent simulation continued from the final state of the partially coupled run. We disregarded the first 30 years to allow for model spin-up and adjustment (chosen based on the ocean energy diagnostic) and analysed the last 100 yr of each simulation. Below we describe the results from the moist model runs only; in 445 the present parameter regime, the qualitative and quantitative behaviour of the companion dry model turned out to be very similar and hence the dry-model results are not shown here (see section 6 for further discussion).
The atmospheric mean state does not appear to depend substantially on whether the temperature feedback on a.m.l. wind is included in the model or not. In each case the atmosphere is characterized by a straight climatological jet with a reasonable vertical shear (Fig. 1) and some zonal modulation; in particular, surface winds tend to be a bit stronger over ocean (due to 450 reduced surface drag), but barotropic wind exhibits an opposite modulation (weaker wind over ocean), consistent with reduced temperature gradient over ocean (Fig. 2). On top of this mean state, the atmosphere is characterized by a vigorous synoptic turbulence (Fig. 3).
The time-mean ocean currents (Fig. 4) represent large-scale subtropical and subpolar gyres and strong inertial recirculations, which help maintain an intense eastward jet in the control and partially coupled simulations. The inertial recirculations 455 largely collapse in the fully coupled simulation (see, for example, the discussion in section 3 and Hogg et al., 2009), with barotropic transports there (~40 Sv) only about 1/3 of those in the control and the partially coupled runs. Accordingly, the | F s ′ | eastward jet becomes much weaker and so is the climatological SST front (Fig. 5); this has probably much to do with anomalous Ekman pumping structures over the inertial gyres seen in the fully coupled simulation (Fig. 5, bottom right). The ocean resolution (of 10 km) requires a relatively high viscosity (and low Reynolds number), and the ocean state can probably 460 be characterized as eddy-permitting, rather than eddy resolving here (Fig. 6). The current simulations are only meant to illustrate a qualitative performance of the model. Simulations at higher resolutions (5 km) and Reynolds numbers may be advisable in all cases (cf. Martin et al. 2020).
We now show the moist characteristics of the model in the control simulation with 'dry entrainment' formulation (23) The evaporation (Fig. 10, left) is prescribed and zonally uniform over land (see section 4) but is active over ocean, exhibiting reduced values to the north and enhancement along the southern and western boundary of the ocean and the double-gyre confluence zone. Atmospheric boundary-layer precipitation is also enhanced over these areas (Fig. 10, right) and exhibits 475 relative minima over the rest of the ocean. Globally, precipitation reaches local maximum at the southern boundary of the model and global minimum at the northern boundary and has a dipolar zonal structure around the axis of the channel, with precipitation minimum/maximum at the anticyclonic(equatorward)/cyclonic(poleward) flanks of the mid-latitude jet, respectively. These features also translate, to some extent, to the precipitation distribution in the atmospheric interior (Fig.   11), although land-sea contrast in precipitation in the interior is opposite in sign to the one within the atmospheric boundary 480 layer. The climatological distribution of precipitation is the result of averaging over intermittent, in space and time, precipitation episodes, as illustrated by a snapshot example in Fig. 12.
Finally, we present here initial evidence for an important effect of the temperature-dependent wind-stress formulation on the low-frequency dynamics of MQ-GCM. This effect can be noticed in the behaviour of the leading EOF of SST (Fig. 13), which is dominated, in all simulations, by a monopolar (in y) SST pattern in the region of ocean's eastward jet and its 485 extension. The intensity, meridional localization and west-to-east scale of this pattern are the largest in the partially coupled run, which has the strongest oceanic turbulence capable of affecting mesoscale winds above the oceanic eastward jet (Fig.   13, top middle panel). This variability tends to be suppressed in both the control run (with no direct SST effect on the a.m.l. winds: Fig. 13, top left panel) and the fully coupled run with active mesoscale coupling (but in which the ocean eddies are partly suppressed: Fig. 13, top right panel), indicating the importance of both the ocean eddies and air-sea mesoscale 490 coupling for this mode. Furthermore, this mode's time-dependence is characterized by a pronounced interdecadal variability in the partially coupled simulation, whereas the dominant time scales in the control and fully coupled simulations are shorter (interannual-to-decadal) and the associated variances -smaller (Figs. 13, bottom row), with the energy-density ratio of the partially coupled to control run of about 6 at low frequencies. It is not immediately clear, however, whether this mode imprints itself onto the atmosphere even in the partially coupled simulation. There are indications that this run's leading jet-495 shifting EOF of the atmospheric streamfunction (analogous to that of the control run shown in the top panel of Fig. 14) has an enhanced energy at the low-frequency end of the spectrum compared to the control and fully coupled simulations (Fig. 14, bottom), but this enhancement is not statistically significant and may be due to sampling. Longer and -most importantly -more highly resolved simulations, in both the ocean and the atmosphere (cf. Martin et al. 2020), are required to gauge the potential of the active mesoscale air-sea coupling to fuel decadal climate modes. 500

Discussion and conclusions
Note that the dynamical core of the present MQ-GCM model is identical to that of the original Q-GCM model, which has already been used in a suite of studies addressing mid-latitude climate variability, while the newly added physics elements have either been tested and verified in ocean-only or atmosphere-only settings (e.g., temperature dependent wind stress: Feliks et al. 2004Feliks et al. , 2007Hogg et al., 2009) or are largely analogous, in their numerical formulation, to the previous Q-GCM 505 elements (e.g., moisture advection and time stepping scheme is analogous to that of the mixed-layer temperature), or are borrowed from similar ocean-only or coupled models (e.g., evaporation and latent heat exchange parameterizations borrow from Kravtsov et al. 2006and Deremble et al. 2012, 2013. The novelty of the present model is in that all these elements are brought together in a fully coupled setting, which makes it a unique numerically efficient tool for exploring possible dynamics of the midlatitude coupled climate variability. While of intermediate complexity, the model is still fairly involved 510 and no reference analytical solutions to directly verify the accuracy of its numerical implementation are available. Furthermore, it is well known from decades of numerical experimentation with analogous ocean-only models (see, e.g., section 1 and references therein) that the regimes of simulated geostrophic turbulence strongly depend on the model's effective Reynolds number, with larger values of this number achievable at a higher horizontal resolution. An analogous statement pertains to frontal air-sea interactions, which may be sensitive to the horizontal resolution in the atmospheric 515 component of the model (see Feliks et al. 2004Feliks et al. , 2007. Future studies of the model sensitivity to the horizontal resolution, in the parameter regimes with high Reynolds numbers and mesoscale-resolving atmosphere are perhaps the ones to bring about the most interesting dynamical insights.
In the present model development effort, we only utilized eddy-permitting ocean resolution and a coarse-grid atmospheric model, which may be behind the similarity between the dry-model and moist-model behaviour noted in section 5. Overall, 520 the MQ-GCM model exhibits rich moisture dynamics reminiscent of many of the observed properties of the Earth's climate system. It remains to be seen whether these and other processes (such as mesoscale air-sea coupling) affect in fundamental ways the dynamics of the mid-latitude ocean-atmosphere system's coupled decadal variability. Preliminary results above indicate that the model's low-frequency variability is indeed sensitive to the details of air-sea interaction; furthermore, both dry and moist versions of the atmospheric model -in parameter ranges corresponding to strong thermal driving and 525 intermediate surface drag (e.g., ) -now exhibit bimodality of the type documented in Kravtsov et al. (2005Kravtsov et al. ( , 2006Kravtsov et al. ( , 2007Kravtsov et al. ( , 2008, which is likely to be important for any decadal climate modes supported by the model (not shown here).

Summary of model updates and code modifications
Changes to the model formulation are summarized below in Table 2. Table 3 outlines the corresponding changes in the Q-530 GCM source code.

Appendix: QGPV equation in a layered atmospheric model
Consider an ideal-gas dry atmosphere comprised of layers with constant potential temperatures . Using the definition of potential temperature, ideal gas law, assuming hydrostatic balance and dropping, in this section, the left superscript "a" that 535 denotes atmospheric variables in the main text -(A1) -one can express the pressures within each layer as The perturbation fields can be found by requiring the pressure to be continuous across each atmospheric interface, 540 is the pressure at the top of the atmospheric mixed layer. For example, from (A1) -(A3) we have, for the first atmospheric layer (A4a) 545 Here, the lower layer streamfunction is given by 550 (A5b) where K is the representative atmospheric potential temperature taken here to be equal, approximately, to the vertical average of individual-layer potential temperatures. In an analogous way, we find, for streamfunctions in layers 2 and 3 (A5c) 555 (A5d) where we assumed that the differences between the potential temperatures of individual layers are small compared to and estimated reduced gravities as (A6) From (A5) it follows that 560 (A7) and that the "dynamic pressures" in Hogg et al. (2014) are equal to .
Furthermore, under synoptic scaling, the continuity equation in each layer is approximately but the last term on the right-hand side of (A8) -approximating -is smaller than other terms and can be neglected 565 under QG scaling, since km and is on the order of the Rossby number.
Thus, the above scaling arguments and calculations demonstrate the validity of Boussinesq approximation for quasigeostrophic compressible atmosphere and justify the uniform treatment of oceanic and atmospheric dynamics in Hogg et al. (2014).

Code availability
The updated code alongside basic instructions on its use (see readme file there), as well as restart files for the six simulations described in this paper, are publicly available from GitHub at https://github.com/GFDANU/q-gcm (Kravtsov et al. 2021).
To use the code, one should replace the routines and scripts of the original source code (publicly available at http://www.q-575 gcm.org with the reference to https://github.com/GFDANU/q-gcm) summarized in Table 3 by their updated versions (contained in the new MQ-GCM folder on GitHub) and compile/run the resulting executable in the same way as before (see Hogg et al., 2014).

Author contribution
SK initiated the study, led the development of MQ-GCM modifications, ran some of the numerical experiments and 580 produced the initial draft of the paper; IM performed and analysed coupled and uncoupled experiments with and without mesoscale air-sea coupling using the original and updated versions of Q-GCM; all co-authors contributed to MQ-GCM development, to analysing the numerical simulations and to writing of the manuscript. We dedicate this paper to the memory of our friend and colleague Prof. Peter Killworth, one of the developers of the original Q-GCM model. Table 1: Mean-state parameters derived from (7)-(11), except for the last two rows detailing the ocean mean state based, loosely, on the observed oceanic vertical structure (note the difference here with the values used in the previous Q-GCM version).