Abstract. For the class of quantum integrable models generated from the q−Onsager algebra, a basis of bispectral multivariable q−orthogonal polynomials is exhibited. In a first part, it is shown that the multivariable Askey-Wilson polynomials with N variables and N + 3 parameters introduced by Gasper and Rahman [1] generate a family of infinite dimensional modules for the q−Onsager algebra, whose fundamental generators are realized in terms of the multivariable q−difference and difference operators proposed by Iliev [2]. Raising and lowering operators extending those of Sahi [3] are also constructed. In a second part, finite dimensional modules are constructed and studied for a certain class of parameters and if the N variables belong to a discrete support. In this case, the bispectral property finds a natural interpretation within the framework of tridiagonal pairs. In a third part, eigenfunctions of the q−Dolan-Grady hierarchy are considered in the polynomial basis. In particular, invariant subspaces are identified for certain conditions generalizing Nepomechie's relations. In a fourth part, the analysis is extended to the special case q = 1. This framework provides a q−hypergeometric formulation of quantum integrable models such as the open XXZ spin chain with generic integrable boundary conditions (q = 1).MSC: 81R50; 81R10; 81U15; 39A70; 33D50; 39A13.