<p>Thermotactic microorganisms move toward warmer regions when a temperature gradient is present, and this behavior can generate complex flow patterns known as bioconvection in fluid-saturated porous media. The main objective of this work is to present the linear and nonlinear (energy) stability analysis of thermotactic bioconvection in a porous medium and to quantify the subcritical instability region arising from finite-amplitude disturbances. The novelty of the study lies in the application of a rigorous energy method to thermotactic microorganism suspensions in porous media, together with a systematic comparison between linear stability thresholds (including both stationary and oscillatory modes) and nonlinear energy stability limits corresponding to stationary disturbances. The linear stability analysis employs the normal-mode approach to determine the critical conditions at which convection first appears. The nonlinear stability analysis is conducted using the energy method, in which an energy functional is constructed and its decay is used to identify the nonlinear stability limit. A single-term Galerkin approximation is employed to solve the resulting variational problem and to estimate the nonlinear critical Rayleigh number. The results show that increasing the swimming speed (<i>Pe</i>) of microorganisms destabilizes the system and causes bioconvection to occur at lower threshold values. However, increasing the porosity (<InlineEquation ID="IEq1"> <EquationSource Format="TEX">\(\epsilon \)</EquationSource> <EquationSource Format="MATHML"><math> <mi>ϵ</mi> </math></EquationSource> </InlineEquation>) of the medium has a stabilizing effect. The thermal Rayleigh number (<i>Ra</i>) enhances convection, while the modified Darcy number (<InlineEquation ID="IEq2"> <EquationSource Format="TEX">\(\tilde{Da}\)</EquationSource> <EquationSource Format="MATHML"><math> <mover accent="true"> <mrow> <mi mathvariant="italic">Da</mi> </mrow> <mo stretchy="false">~</mo> </mover> </math></EquationSource> </InlineEquation>) exhibits a dual role: it stabilizes the system at higher wavenumbers but destabilizes it at lower wavenumbers. A distinct subcritical region of instability is observed, and this region shrinks as microorganism motility increases. Also, oscillatory bioconvection is investigated using linear stability analysis. It is found that higher microorganism motility significantly reduces the oscillatory instability threshold, whereas increasing porosity, stable thermal stratification, and larger values of <InlineEquation ID="IEq3"> <EquationSource Format="TEX">\(\tilde{Da}\)</EquationSource> <EquationSource Format="MATHML"><math> <mover accent="true"> <mrow> <mi mathvariant="italic">Da</mi> </mrow> <mo stretchy="false">~</mo> </mover> </math></EquationSource> </InlineEquation> delay the onset of time-periodic convection and stabilize the system. These results are useful for applications such as porous bioreactors, tissue engineering scaffolds, and systems designed to control microbial transport.</p>