This article reports on comprehensive assessment regarding conjoined effect of discretization order of the continuous Boltzmann equation and discrete forcing scenario in the feasibility of lattice Boltzmann method (LBM) as an interlinked hydrothermodynamics solver. A two-dimensional natural convection phenomena within the geometrical pore arrangement of the concrete structure was modeled as a synergetic embodiment of the discrete lattice fluid and thermal counterparts. Four different combination strategies of LBM implementation were carefully examined, elucidating theoretical and numerical segments, in order to expose any plausible discrepancy in the retrieved steady-state solutions. It was found that the numerical outcomes from distinct LBM strategies return equivalent results upon few selected key physical properties, despite incisive difference that emerges in the theoretical aspect. Excellent agreement with classical computational techniques was observed for all considered treatment options, underlining validity of our strategies. This study represents a further step toward clarification of the convoluted issue regarding proper selection of discretization order and forcing scheme in LBM simulation.