We present an efficient computational approach for simulating component transport within single-phase free flow in the boundary layer over porous media. A numerical model based on this approach is validated using experimental data generated in a climate-controlled wind tunnel coupled with a soil test bed. The developed modeling approach is based on a combination of the lattice Boltzmann method (LBM) for simulating the fluid flow and the mixed-hybrid finite element method (MHFEM) for solving constituent transport. Both those methods individually, as well as when coupled, are implemented entirely on a GPU accelerator in order to utilize its computational power and avoid the hardware limitations caused by slow communication between the GPU and CPU over the PCI-E bus. In order to utilize vast computational resources available on modern supercomputers, the implementation is extended for distributed multi-GPU computations based on domain decomposition and the Message Passing Interface (MPI). We describe the mathematical details behind the computational method, focusing primarily on the coupling mechanisms. The performance of the solver is demonstrated on a modern high-performance computing system. Flow and transport simulation results are validated and compared herein with experimental velocity and relative humidity measurements made above a flat partially saturated soil layer exposed to steady air flow. Model robustness and flexibility is demonstrated by introducing cuboidal bluff-bodies to the flow in several different experimental scenarios. The experimentally measured values are available in a publicly available dataset that can serve as a benchmark for future studies. Finally, we discuss potential improvements for the model as well as future experimental efforts.