We perform finite element modeling of pulse-echo ultrasound of a hard inclusion in a soft tissue to gain a better understanding of B-mode image brightness characteristics. We simulate a pressure wave emitted by an ultrasound transducer through the inclusion-tissue medium by prescribing suitable boundary conditions, and collect the scattered wave response to simulate the behavior of the transducer array used for pulse-echo ultrasound. We form B-mode images from simulated channel data using standard delay and sum beamforming. We establish the accuracy of the finite element model by comparing the point spread function with that obtained from Field II ultrasound simulation program. We also demonstrate qualitative validation by comparing the brightness characteristics of rough and smooth surfaced circular inclusions with experimental images of a cylindrical metal tool immersed in a water tank. We next conduct simulation studies to evaluate changes in B-mode image brightness intensity and contrast related to different constitutive properties, namely, compressibility of the inclusion, impedance contrast between the host and inclusion, and surface roughness of the inclusion. We find that the intensity observed behind a hard inclusion in the axial direction is strongly affected by the compressibility and roughness of the inclusion. Also, the perceived width of the stone based on the intensity is greater for rougher stones. Our study indicates that imaging of compressible inclusions may benefit from targeted B-mode image forming algorithms. Our modeling framework can potentially be useful in differentiating hard inclusions from surrounding parenchyma, and for classifying kidney stones or gallstones.