Solving wave equations is pivotal for various downstream applications such as microseismic imaging and CO2 monitoring. However, obtaining such solutions is costly and prone to errors, requiring a fast and flexible alternative simulation tool to address the limitations. Knowledge-guided neural networks have emerged as promising surrogate modeling tools. But it is crucial to integrate geophysical domain knowledge to realize their full potential in seismic modeling and subsurface monitoring. This dissertation focuses on developing neural network based PDE solvers (neural PDE solvers) for this purpose. The research begins with physics-informed neural networks (PINNs), which offer several advantages, including adaptability to irregular geometries, meshless solutions, and ease in handling complex coupled equations. However, we found that due to the low-frequency bias and the complex loss landscape of PINNs training, many challenges exist, e.g., the convergence of training, poor representations for high-frequency and multi-frequency wavefields. To tackle these challenges, I developed a modified PINN incorporating positional encoding to mitigate the low-frequency bias issue. To further make the training for PINN faster, I developed an even more efficient PINN using hash encoding, which offers a locally aware (at multi-resolution) coordinate input to the neural network (NN). Recognizing the importance of network architecture, I also incorporated prior knowledge of the seismic wavefield into the network design using Gabor basis function-incorporated neural networks. To improve the convergence and accuracy for high-frequency wavefields, I developed a PINN evolution based on frequency upscaling and neuron splitting (PINNup), and introduced a single reference frequency loss function to enhance multi-frequency wavefield representation. Using these advancements in PINNs, I formulated a direct microseismic source imaging framework with hard constraints. To support instant simulation across various velocities, I also developed a learned frequency-domain scattered wavefield solution using a neural operator. Applying this neural PDE solver directly to the full waveform inversion (FWI) will result in noisy inverted velocities. Thus, I further enhanced the approach with a physics-informed waveform inversion to tackle this issue. To monitor CO2 injections and predict the CO2 plume evolution, I utilized a video-diffusion-based framework for end-to-end multiphysics monitoring and forecasting, considering the complexity of real-world subsurface monitoring.