Thermal convection in a planar domain with periodic horizontal extent bounded below by a heated rigid plate and above by a deformable free surface with temperature dependent surface tension is numerically simulated under the effect of gravity in a variational formulation. In the horizontal direction, the assumption of periodicity allows the use of Fourier series expansion. In the vertical direction, Legendre Lagrangian interpolants expansion are used for the velocity and temperature, while orthogonal Legendre polynomial expansion with order less than that for velocity is used for the pressure [1]. A lower order expansion for pressure eliminates the problem of pressure boundary condition at the lower plate. Since the governing equations including the free boundary are nonlinear, a Newton-like iterative method is used based on successive total linearization with respect to the velocity, temperature, pressure and the position of the free boundary. The linearization is performed on the continuous form of the variational integrals over the domain $\Omega(h)$ with the free boundary $\Gamma_f(h)$ as shown below for a flow quantity $q$: $$\int_{\Omega(h)}^{}q\,d\Omega \approx \int_{\Omega(h^{(n)})}(q^{(n)}+\delta q)\,d\Omega + \int_{\Gamma_f(h^{(n)})}{\delta h}\,q^{(n)}\,ds$$ at each iteration step $(n)$ so that $q^{(n+1)}=q^{(n)}+\delta q$ and $h^{(n+1)}=h^{(n)}+\delta h$ are updated where $\delta h$ is the variation of the free-surface along the normal $\textbf{n}$ [2]. High order polynomial expansions and the Newton-like iteration procedure lead to accurate results and fast-converging iterations. REFERENCES [1] M. R. Schumack, W. W. Schultz, J. P. Boyd. Spectral method solution of the Stokes equations on nonstaggered grids", J. Comp. Phys., Vol. 94, pp. 30-58, 1991. [2] N. P. Kruyt, C. Cuvelier, A. Segal, J. van der Zanden, Total linearization method for solving viscous free boundary flow problems by the Finite Element Method. Int. J. Num. Meth. Fluids, Vol. 8, pp. 351-363, 1988.