The authors consider an incompressible fluid flowing in the 2D domain \( [0,\pi A]\times \lbrack 0,\pi ]\) and which is heated at a temperature higher at the bottom boundary than at the top boundary.Considering the Boussinesq approximation, they write the equations \(\partial _{t}u+u\cdot \nabla u=-\nabla p+\sigma \nabla ^{2}u+\sigma RaT\widehat{z}\), \(\nabla \cdot u=0\), \( \partial _{t}T+u\cdot \nabla T=\nabla ^{2}T\), where \(u\) is the velocity, \(p\) the pressure, \(T\) the temperature, and \(Ra\) Rayleigh number. Introducing the stress function \(\psi \) such that \(u=(\partial _{z}\psi ,-\partial _{x}\psi ) \) they end with the dimensionless Boussinesq equations \(\partial _{t}\nabla ^{2}\psi -\{\psi ,\nabla ^{2}\psi \}=\sigma \nabla ^{4}\psi +\sigma \partial _{x}\theta \), \(\partial _{t}\theta -\{\psi ,\theta \}=\nabla ^{2}\theta + \mathcal{R}\partial _{x}\psi \), where \(\{f,g\}=\partial _{x}f\partial _{z}g-\partial _{z}f\partial _{x}g\), and \(\theta =\pi \mathcal{R}(T_{c}-T)\), with \(T_{c}=1-z/\pi \). The boundary conditions \(\psi ,\partial _{z}^{2}\psi ,\theta =0\) at \(z=0,\pi \) are imposed. The authors compute a value of the Nusselt number \(Nu\) in terms of infinite-time averages of \(\psi \) and \( \theta \). The next step consists to introduce projections of the above Boussinesq equations on Fourier modes leading to the expressions \(\psi (x,z,t)=\psi _{11}(t)\sin (kx)\sin (z)+\psi _{12}(t)\cos (kx)\sin (2z)+\psi _{01}(t)\sin (z)+\psi _{03}(t)\sin (3z)\), \(\theta (x,z,t)=\theta _{11}(t)\cos (kx)\sin (z)+\theta _{12}(t)\sin (kx)\sin (2z)+\theta _{02}(t)\sin (2z)+\theta _{04}(t)\sin (4z)\), with \(k=2/A\), then to a coupled system of ODEs for the coefficients of these expressions. The authors compute three steady state solutions of this system corresponding to Rayleigh numbers \(\mathcal{R}_{L_{1}}=(k^{2}+1)^{3}/k^{2}\), \(\mathcal{R} _{L_{2}}=(k^{2}+4)^{3}/k^{2}\) and \(\mathcal{R}_{TC_{1}}\) given in terms of \( \mathcal{R}_{L_{1}},\sigma ,k\), with bifurcation structures that they illustrate with figures. They use MATLAB to compute time-dependent states. The main purpose of the paper is to build upper bounds for the optimal time averaged heat transport among all trajectories of the system of ODEs. The authors use a general method for ODEs written as \(\overset{.}{x}=f(x)\) where \(f:\mathbb{R}^{n}\rightarrow \mathbb{R}\) is continuously differentiable and a quantity \(\phi :\mathbb{R}^{n}\rightarrow \mathbb{R}\) which is continuous. They define \(\overline{\phi }^{\ast }=\sup_{x_{0}\in \mathbb{R} ^{n}}\overline{\phi }(x_{0})\), where \(\overline{\phi }(x_{0})\) is the infinite-time average of the trajectory \(x(t)\) which starts from \(x_{0}\). They introduce an auxiliary function \(V:\mathbb{R}^{n}\rightarrow \mathbb{R}\) which is continuously differentiable and such that \(\overline{f\cdot \nabla V }=0\) on every trajectory and they obtain the bound \(\overline{\phi }^{\ast }\leq \inf_{V\in C^{1}}\sup_{x\in \mathbb{R}^{n}}[\phi (x)+f(x)\cdot \nabla V(x)]\). Finally, choosing \(\phi =1+\frac{1}{\mathcal{R}}(2\theta _{02}+4\theta _{04})\) so that \(\overline{\phi }=N\), the Nusselt number, they build numerical upper bounds and they analyze the dependence on different parameters. They also compute analytical upper bounds introducing appropriate quadratic auxiliary functions and they illustrate these results with figures.