Interactive comment on “ A 24-variable low-order coupled ocean – atmosphere model : OA-QG-WS v 2 ” by S . Vannitsem

1. You are right. More information should be added concerning the results obtained with previous low-order coupled ocean models. This further highlights the similarities between the present work and the previous ones, with a marked difference with the results of Nese and Dutton (1993) and much more similarities with the work of van Veen (2003). In the former, the activation of the dynamics within the ocean leads to an increase of predictability. This feature contrasts with our results but could probably be

By definition, these low-order models are built in such a way to simplify as far as possible the system under investigation and keep only the key ingredients of interest, as for instance the analysis of the impact of an orography on the instability of atmospheric flows as in Charney and Straus (1980).When dealing with climate the same procedure can be performed by focusing on one specific aspect, for instance the global energy balance of the earth assuming that the dynamics at smaller space and timescales could be modeled based on stochastic processes (Nicolis and Nicolis, 1979).When one is interested in keeping key ingredients of processes acting at very different scales, the problem becomes Published by Copernicus Publications on behalf of the European Geosciences Union.more involved and only a few such models have been developed.A popular approach consists in coupling two low-order models and modifying artificially the typical timescale of one of them (e.g., Goswami et al., 1993;Pena and Kalnay, 2004).This approach could indeed provide an easy way to build such multiscale models, but one loses physical significance.Another interesting model built in this form was proposed by Roebber (1995), in which the low-order Lorenz (1984a)' model is coupled with an oceanic three-box model (with six ordinary differential equations for temperature and salinity) developed by Birchfield (1989), using empirical relations for heat fluxes.This led to a coupled model of nine prognostic variables, with two specific timescales, one for the atmosphere and the other for the ocean.
The other approach consists in starting from a detailed coupled model and systematically reducing the number of modes of the different components.A first attempt made by Lorenz (1984b) led to a coupled ocean-atmosphere loworder model incorporating many processes like condensation, evaporation, and radiative transfer.However, the ocean was only considered as a heat bath.This model was subsequently modified by Nese and Dutton (1993), in which oceanic transport is incorporated in a way similar to Veronis (1963).The final version of this model contains 31 prognostic variables and several diagnostic relations.The coupled model developed by Nese and Dutton (1993) was used to evaluate the impact of the ocean transport on the predictability of the coupled system.They have found that when the ocean dynamics is activated, the predictability as measured by the Lyapunov exponents is increased.Another interesting model developed by van Veen (2003) and derived from first principles combines the three-variable atmospheric system of Lorenz (1984a) and the four-variable ocean model of Maas (1994).In this seven-variable model, a clear distinction between three different timescales is made, one for the atmosphere, one for the deep ocean and one for the ocean surface layer.In this model, a systematic bifurcation analysis has been undertaken and compared with the bifurcation structure of the atmospheric component only.In particular it was shown that the ocean plays an important role close to the bifurcation points of the model, but much less in the chaotic regime.In the latter case the ocean integrates the rapid fluctuations of the atmosphere in a quite passive manner without providing a strong feedback toward the atmosphere.In addition, only single oceanic gyres can develop.
Building on the latter stream of ideas, Vannitsem (2014) proposed to couple two low-order models for the atmosphere and the ocean, derived from quasi-geostrophic equations.This model is intermediate between the "very loworder" coupled models proposed by van Veen (2003), and the more sophisticated process-oriented low-order coupled models of Lorenz (1984b) and Nese and Dutton (1993).It is based on the low-order quasi-geostrophic model of Charney and Straus (1980) and the shallow water quasigeostrophic model of Pierini (2011).The latter is able to simulate the dynamics of single or double oceanic gyres, typical in the Northern Atlantic and Pacific.The coupling is done through momentum transfer at the interface, only.This model has the advantage to be derived from first principles as in van Veen (2003) and Lorenz (1984b), but focusing only on the coupled dynamics associated with the momentum forcing between the two components.It will be referred to as OA-QG-WS v1 (OA-QG-WS for Ocean-Atmosphere Quasi-Geostrophic Wind Stress).An extension has also been proposed in Vannitsem (2014), by adding atmospheric modes as in Reinhold and Pierrehumbert (1982).This second version of the model, whose dynamics was only slightly touched upon in Vannitsem (2014), is the central subject of the present paper, and will be referred to as OA-QG-WS v2.
The degree of sophistication of this low-order model is such that it is not straightforward to evaluate all the coupling coefficients (and the coefficients of the oceanic part), due to the presence of different orthogonal basis functions and inner products for both climate components.These are therefore made available here and some validation test cases are provided for subsequent use of the model by the atmospheric and climate communities.The revision of the model also allowed correcting a few coefficients of the first model version presented in Vannitsem (2014), without qualitative modifications of the results and conclusions.In addition, a few results concerning the dynamical instability of the system are provided, and similarities and dissimilarities with the trends already found in Vannitsem (2014) are discussed.
The original partial differential equations of the model and the choice of the orthogonal modes are presented in Sect. 2. Section 3 is devoted to some properties of the model that could serve as a benchmark.The appendix contains all the coefficients of the model, as described in Sect. 2. In Sect.4, some conclusions are drawn.

The atmospheric model
The atmospheric model, developed by Charney and Straus (1980) and subsequently extended by Reinhold and Pierrehumbert (1982), is a two-layer quasi-geostrophic flow defined on a beta plane.The equations in pressure coordinates are where ψ 1 , ψ 3 , ω are the streamfunction fields at 250 and 750 hPa, and the vertical velocity (i.e., dp/dt), respectively.f 0 is the Coriolis parameter at latitude φ 0 , β = df/dy at φ 0 , σ = −R/p ∂T ∂p − RT pc p is the static stability (where T is the temperature, R the gas constant and c p the heat capacity at constant pressure) considered as constant.k d and k d are the coefficients multiplying the surface friction term and the internal friction between the layers, respectively.(ψ 1 − ψ 3 ) * is a constant thermal forcing of the atmosphere (Newtonian heating).An additional term has been introduced in this system in order to account for the presence of a surface boundary velocity of the oceanic flow defined by (see next section).This would correspond to the Ekman pumping on a moving surface and is the mechanical contribution of the interaction between the ocean and the atmosphere (e.g., Deremble et al., 2012).
Note also that the heating term has not been modified even if heating is coming mostly from the ocean.It is assumed that this heating is a fast process as compared to the dynamics of heat transport in the ocean, thereby transferring almost instantaneously the energy toward the atmosphere.This strong assumption allows isolating the impact of wind-driven interactions between the ocean and the atmosphere.This assumption could be relaxed in a future version of the model in a similar way as in van Veen (2003) or Deremble et al. (2012).
These equations are then adimensionalized by scaling x = x/L and y = y/L, t by f −1 0 , ω by f 0 p and ψ by L 2 f 0 and the parameters are then also rescaled as The fields are expanded in Fourier series over the domain y = [0, π ] and x = [0, 2π/n], and only 10 modes, F k , are retained, obeying the boundary conditions ∂F k /(∂x ) = 0 at y = 0, π .n is the aspect ratio between the lengths of the domain in y and in x, n = 2L y /L x = 2πL/(2πL/n).These modes are and the fields are then expressed as where θ = (ψ 1 − ψ 3 )/2 and ψ = (ψ 1 + ψ 3 )/2.Using the usual inner product, one gets the set of equations reported in the Appendix of the paper of Reinhold and Pierrehumbert (1982) and in Reinhold and Pierrehumbert (1985), leading to 20 ordinary differential equations for the dependent variables ψ k and θ k .The dynamics of this atmospheric model has also been explored with emphasis on the predictability of the atmosphere in the presence of weather regimes in Trevisan et al. (2001).
The presence of the ocean is felt through the coupling associated with the motion of the ocean surface, k d ∇ 2 where is the streamfunction of the oceanic flow as defined in the next section.It is also projected on the different atmospheric modes using the inner product of Eq. ( 4).The coefficients are given in Appendix B.
Note that the thermal forcing term is fixed as in Charney and Straus (1980) and Reinhold and Pierrehumbert (1982) in which the only nonzero term is θ * 1 , which will be referred to as θ * in the sequel.This corresponds to a thermal forcing only dependent on the latitude with a larger contribution in the southern part of the domain.

Ocean model
The ocean model is based on the reduced-gravity quasigeostrophic shallow water model (Vallis, 2006).The basic assumptions behind this equation are that (i) the ocean dynamics can be described by a shallow water fluid layer superimposed over a quiescent deep fluid layer, (ii) the Rossby number Ro = U/(f 0 L) is small, and (iii) the space scale of the process under investigation should not be significantly larger than the deformation radius (typically of a few hundred kilometers for a fluid layer depth of the order of 100 m).The forcing is provided by the wind generated by the atmospheric component of the coupled system.The equation reads where is the velocity streamfunction (or pressure), ρ the density of water, h the depth of the fluid layer, L R the reduced Rossby deformation radius, r a friction coefficient at the bottom of the fluid layer, and curl z τ the vertical component of the curl of the wind stress.Usually in low-order oceanic modeling the latter is provided as an ideal profilein the meridional direction (e.g., Simonnet and Dijkstra, 2002).In the present work, this is provided as a "real" wind field generated by the atmospheric low-order model.Assuming that the wind stress is given by where u and v are the horizontal components of the lower layer geostrophic wind, −∂ψ 3 /∂y and ∂ψ 3 /∂x, respectively, and U and V the corresponding quantities in the ocean, one gets Here the wind stress is proportional to the relative velocity between the flow in the ocean layer and the wind.This slight modification as compared with the version model OA-QG-WS v1 in which the stress was only based on the absolute wind velocity, has been made in order to avoid spurious forcings when the velocities in the atmosphere and the ocean are similar.It is however a correction which is quite marginal in view of the (typically) small amplitudes of the flow field in the ocean.
Using the same domain and the same nondimensionalization procedure as in the atmospheric model, one gets where Let us now define the truncated basis functions on which the streamfunction field is projected.Several truncations were proposed in the literature from two-mode (Jiang et al., 1995) up to four-mode truncations (Simonnet et al., 2005;Pierini, 2011), the latter approach allowing for chaotic behaviors.In the present work, we use the following set of modes, in order to get the free-slip boundary conditions (and no normal flow to the wall) in the domain over which the flow is defined at x = 0, 2π/n and y = 0, π .In addition a specific inner product is adopted for the oceanic model in a similar way as in Pierini (2011), Introducing the truncated fields, m A m φ m , for m = 1, 4, into Eq.( 7) and projecting on each mode using the inner product Eq.( 9), one gets a set of four ordinary differential equations for the variables A m , whose coefficients are all provided in Appendix A. The forcing tendencies, f (m), m = 1, 4, associated with the wind stress as defined by Eq. ( 6), are given by whose coefficients are provided in Appendix B, where should not be confused with the Coriolis parameter f 0 and the parameter f 1 of Appendix A.

Estimation of the main parameters
The estimation of the main physical parameters is made as follows.For the atmosphere, the parameter k is related to the surface drag felt by the lower layer of the two-layer QG model.This is estimated based on the Ekman layer theory (p. 115, Vallis, 2006) as after dividing by f 0 , and where D and d are the thickness of the lower atmospheric layer and the thickness of the Ekman surface layer, respectively.Typically D is of the order of 5000 m and d of the order of 100-1000 m.This implies that k falls in a range of [0.01, 0.1].Here the value is fixed to k = 0.02 (and the other dissipation parameters are fixed to h = k = 2k).For parameter δ, one can use the estimate done by Nese and Dutton (1993).The dimensional forcing coefficient is given by where ρ a and ρ o are the densities of the air and of the sea water, respectively.h is the thickness of the ocean layer and C D the surface friction coefficient.With C D ≈ 0.001,  6) is equivalent to For the thermal forcing, the same approach as in Charney and Straus (1980) and in Reinhold and Pierrehumbert (1982) is adopted, through the use of the thermal wind relation.θ * is therefore allowed to vary from [0, 0.2].

Results of the integration of OA-QG-WS v2
In this section, some statistical and dynamical properties of the model are reported as a benchmark.The numerical scheme used is a second-order temporal scheme known as the Heun scheme (see Kalnay, 2003) with a time step of 0.01 time unit.The parameter values used are listed in Table 1, while the behavior of the system is explored by varying δ and θ * .The dimensional time unit is equal to 0.11215 days.

Model trajectories and mean fields
Figure 1 displays the temporal evolution of the A i variables of the ocean component for about 10 years starting after 200 000 days of integration.Interestingly a long-range variability emerges as in Vannitsem (2014).
As already alluded in Vannitsem (2014), this new version of the model allows for the development of double gyres.Figure 2 displays the mean streamfunction fields for different values of the key parameters θ * = 0.077, θ * = 0.10, and θ * = 0.14, after a long integration of about 3.5 × 10 8 days.Two different initial states in phase space are used for θ * = 0.077 in panels a and b.Depending on parameter (and maybe initial state in phase space) choice, different mean configurations and sizes of gyres could develop in the basin.But as reflected in Fig. 1, a large variability on a wide range of timescales is also present around these mean fields leading to a variable transport in the ocean basin.The temporal Table 1.Dimensional and nondimensional parameters used in the coupled ocean-atmosphere model.

Dimensional parameters
Nondimensional parameters variation of these mean values are illustrated in Fig. 3, for θ * = 0.077 and θ * = 0.14, starting from two different initial conditions.The convergence is very slow due to the natural long-term variability of the ocean embedded in this system.The presence of different attractors cannot be confirmed or excluded at this stage, due to the blurring of the large natural variability of the system.This analysis would need even longer model integrations, with a higher-order numerical scheme in order to better control the numerical error as suggested by the anonymous referee.Two codes (in Fortran and Lua) used to integrate the model (with the secondorder Heun method) and compute these averaged quantities are provided as Supplement and can be used freely, provided proper reference to the source is made.
Figure 4 displays the power spectra of modes ψ 1 and A 1 , as obtained using a time series of about 73 500 days for θ * = 0.14 (sampled every 0.56075 days, one point every 500 adimensionalized time steps).The atmospheric field displays a flat spectrum for small frequencies and decays at the large ones.The typical timescale of transition between these two regimes is of the order of 30 days for this large-scale atmospheric mode.For the oceanic mode, the power spectrum is continuously decaying closely following a power law, indicating long-range time dependences (in agreement with the visual inspection of Fig. 1).A change of slope is also visible in this log-log plot, around a timescale of 30 days, reflecting the change of statistical properties in the atmosphere.For low frequencies (between ω = 0.0001 and ω = 0.2, the slope of the decay is close to −2, suggesting a dynamics close to a red noise.For large frequencies, the slope is much sharper with a value close to −4.At low frequencies the ocean acts as an integrator of the "white" noise produced by the atmosphere, by analogy with a Brownian motion or an Ornstein-Uhlenbeck process.

Chaotic dynamics
Sensitivity to initial conditions is one of the main properties of the atmosphere.In dynamical systems theory, this property is usually quantified by evaluating the Lyapunov exponents.These quantities also allow for distinguishing between the typical solutions generated by the system of ordinary differential equations for some specific parameters.For a detailed discussion of these typical solutions and the numerical algorithms used for their evaluation, see Parker and Chua (1989).In short, these quantities characterize the amplification of small amplitude initial condition errors in time and are evaluated in the so-called tangent space of the model trajectory (Legras and Vautard, 1996), formally characterized by the Jacobian matrix of the flow.In this tangent space, it can be shown that there exists a set of (characteristic) vectors, u i (t), i = 1, . . ., n, and a corresponding set of (characteristic) numbers, σ i , quantifying the degree of amplification of small perturbations, δx i (t), along these vectors.These characteristic numbers are known as the Lyapunov exponents and are given by If one of these exponents is positive, then the system is sensitive to initial conditions and the solution is chaotic.If the largest one is 0 and the others negative, then the solution is periodic.If the two largest exponents are 0 and the others negative, the solution lives on a 2-torus.Practically it is not necessary to know these specific vectors, u i (t), i = 1, . . ., n, to get the Lyapunov exponents and any basis of independent vectors can be used, because the amplification of any L-dimensional volume in phase space will amplify on average with a rate equal to the sum of the first L Lyapunov exponents (e.g., Legras and Vautard, 1996).Numerically one uses a basis which is regularly orthonormalized in order to avoid the collapse of all the vectors along the dominant instability direction (e.g., Parker and Chua , 1989).
One of the main properties of this new version of the model is the possibility of having a "large" number of positive Lyapunov exponents, and hence a "large" attractor dimension.Figure 5a displays the variations of the first, second and third Lyapunov exponents as a function of θ * for δ = 0.001938.For values of θ * smaller than 0.055, stable steady states are found with a set of four negative Lyapunov exponents of very small amplitude (e.g., for θ * = 0.02, σ 1 = −0.00128,σ 2 = −0.00128,σ 3 = −0.00133,and σ 4 = −0.00133day −1 ) and the next ones with an amplitude 1000 times larger.At θ * = 0.055, a periodic solution emerges with a first exponent equal to σ 1 = −1.110−8 day −1 .For larger values up to θ * = 0.065, quasi-periodic solutions (2-torus) appear, as well as for parameter values between 0.087 and 0.095.Between 0.065 and 0.087, chaotic solutions separated by small periodic windows are prevailing.Beyond 0.095, the dynamics become chaotic and no periodic solutions were found for the parameter values explored.For large values of θ * the dynamics becomes wilder with a dominant exponent close to σ 1 = 0.50 day −1 for θ * = 0.16, a value larger than the ones found for more realistic synoptic-scale dynamics (Vannitsem and Nicolis, 1997;Snyder and Hamill, 2003).Figure 5b  model has therefore more flexibility since one can easily get different configurations in terms of dynamical instability, by changing the main parameter θ * .A detailed analysis of the transitions from quasi-periodic motions to chaotic behaviors will be investigated in the future as in recent works (Broer et al., 2011;Sterk et al., 2010, among others).
Figure 6 displays the dependence of the amplitudes of the Lyapunov exponents and the number of positive exponents as a function of the coupling parameter δ, for three different values of θ * .As in Vannitsem (2014), the trends of the Lyapunov properties as a function of δ can be very different for different values of θ * .The values of the exponents for θ * = 0.0825 are very sensitive to δ, with sharp transition from (quasi-)periodic solutions to chaotic behaviors around δ = 0.009.This interesting feature suggests that δ plays a crucial role in setting up the transition from nonchaotic to chaotic regimes in the coupled system.A full understanding of this transition should be obtained through a systematic analysis of the bifurcation diagram of this system (and it will be the subject of a future investigation).For θ * = 0.10 and θ * = 0.14 an increase is found for the first two exponents (but very weak for θ * = 0.14), while a third positive one emerges when δ is increased.Interestingly, these results confirm the tendency already reported in van Veen (2003), indicating that the presence of the ocean has a stronger influence on the dynamics of the atmosphere close to the periodic windows.
The sensitivity to δ is also illustrated in Fig. 6d in which the Kolmogorov-Sinai entropy is shown, displaying a systematic increase for the three values explored.These trends are opposite to those discovered in Nese and Dutton (1993).Their results are most probably associated with the way the heat is transported in the ocean basin and then transferred toward the atmosphere in their model, a feature not present in our model.This is worth investigating further in the future by adding thermal exchanges between the atmosphere and the ocean.
For all the cases explored, the number of positive Lyapunov exponents also has a tendency to increase with the amplitude of the coupling δ.This feature is similar to what was found in OA-QG-WS v1, further reflecting the importance of the coupling between the ocean and the atmosphere.
To further understand this increase of instability as a function of the coupling parameter, the mean absolute amplitude of the (backward) Lyapunov vectors along the different variables of the coupled system has been computed.Figure 7 displays the results for the first (backward) Lyapunov vector (see Legras and Vautard, 1996) corresponding to the dominant Lyapunov exponent for the same parameter as in Fig. 5c and for three different values of δ.The first 10 points correspond to the barotropic variables of the system, the next 10 points to the baroclinic ones, and the last 4 points to the ocean variables.Clearly the projections along the atmospheric variables do not change as a function of the coupling δ, contrary to the projection along the ocean variables.A similar picture is found for the other backward Lyapunov vectors.This suggests that the increase of instability is mainly associated with an increase of the projection of the vectors along the ocean variables, and not the baroclinic or barotropic instability within the atmosphere.This conjecture is worth investigating further in the future through a detailed analysis of the bifurcation diagram and of the characteristic vectors (also called covariant vectors) of the system, which are (nonorthogonal) intrinsic directions of instabilities (see Legras and Vautard, 1996).

Conclusions
In this paper, a new version (OA-QG-WS v2) of a low-order coupled ocean-atmosphere model is presented, containing 24 ordinary differential equations.This model describes the dynamics of the large-scale flows at midlatitudes of a baroclinic atmosphere interacting with an ocean layer under wind forcings (or momentum exchanges).This coupled model displays features with strong resemblance with the dynamics found at midlatitudes, with a chaotic dynamics of the atmosphere at short timescales of the order of a day and a decadal variability of the ocean layer.In contrast with the model version OA-QG-WS v1 (Vannitsem, 2014)  ean absolute amplitude of the first (backward) Lyapunov vector along the variables of the om 1 to 10, ψ i , from 11 to 20, θ i , and from 21 to 24, A i , as obtained after an integration of 10 6 37 Fig. 7. Mean absolute amplitude of the first (backward) Lyapunov vector along the variables of the system, from 1 to 10, ψ i , from 11 to 20, θ i , and from 21 to 24, A i , as obtained after an integration of 10 6 days.The other parameters as in Fig. 6c.
attractors (associated with a larger number of positive Lyapunov exponents) can be found, and double gyres can develop in the ocean basin in the presence of a chaotic atmosphere.
The Lyapunov instability properties of the flow have also been explored.Interestingly, for the set of parameters chosen, a transition from periodic to chaotic regimes occurs at a value of the bifurcation parameter close to θ * = 0.065.Close to this value, the dynamics is also highly sensitive to the values of the coupling parameter δ, with a possibility of a sharp transition from periodic to chaotic regimes.For large values of θ * , the dominant exponent is less sensitive to δ, in contrast to the lowest amplitude positive exponent.In addition, the number of positive Lyapunov exponents has a tendency to increase with δ regardless of what θ * is, suggesting an increase of the dimension of its attractor in phase space.The latter characteristic was also found in the first version (OA-QG-WS v1) of the model.
As suggested by the analyses reported above, this new model version is an interesting candidate for subsequent analyses of the dynamical properties of coupled systems.In addition, it can be used for testing tools developed for coupled ocean-atmosphere systems in the context of data assimilation, post-processing, and ensemble forecasting, among others.All the coefficients of the (ocean) model and of the coupling terms are also provided, allowing for an easy implementation.Two codes (in Fortran and Lua) combining the atmospheric and oceanic components are also provided as Supplement.

Coefficients of the ocean component of the model
where

Appendix B
Coefficients of the coupling between the ocean and the atmosphere (1 − e α2π/n ), and where the B i = ψ 3 i = ψ i −θ i are the atmospheric streamfunction variables (mode amplitudes) in the lower layer.

Fig. 2 .
Fig. 2. Average streamfunction field of the ocean for δ = 0.001938 and θ * = 0.077 (a), 0.077 (b), 0.10 (c) and 0.14 (d), as obtained from a long integration of about 3.5 × 10 8 days.Note that (a) and (b) are obtained with the same parameters but different initial states in phase space.

Fig. 6 .Fig. 6 .
Fig. 6.Values of the 3 first Lyapunov exponents as a function of the coupling parameter δ for θ * = 0.0825 (a), 0.10 (b), and 0.14 (c).(d) Variation of the Kolmogorov-Sinai entropy as a function of δ for the same values of θ * .